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Multi-Physics  Modeling  and  Simulation  of  Process-Induced  Stresses  in  Polymer-Matrix 
Composites  /  Dr.  Kishore  Pochiraju 

Summary: 

A  three-dimensional  thermo-chemical  simulation  of  the  curing  phase  of  Liquid  Composite  Molding 
Process  coupled  with  process-induced  stress  and  deformation  analyses  are  developed.  The  mechanics 
model  considers  orthotropic  and  thermo-chemically  varying  viscoelastic  stiffhess  of  the  composite 
material  and  varying  boundary  conditions  during  the  process  cycle.  The  residual  stress  and  deformation 
profiles  resulting  from  several  processing  histories  including  the  effect  of  thermal  expansion,  and 
chemical  shrinkage  are  studied.  Several  numerical  simulations  are  performed  to  compare  and  validate 
the  developed  method  with  the  existing  solutions. 

What  Was  Accomplished: 

A  Galerkin  finite  element  with  an  Jacobi  Conjugate  Gradient  (JCG)  iterative  solver  is  used  to  solve  the 
thermal  equilibrium  and  chemical  kinetics  during  the  curing  phase  of  RTM  process.  The  mold 
temperature  history  and  molding  process  parameters  are  translated  into  initial  and  boundary  conditions 
of  the  simulation.  The  temperature  and  degree  of  cure  fields  and  their  gradients  are  obtained  during  the 
curing  process.  Two  material  constitutive  models  describing  the  thermoelastic  and  viscoelastic  behavior 
of  the  material  are  considered.  Effective  composite  mechanical  properties  were  computed  using  the 
instantaneous  resin  and  fiber  properties  in  a  self-consistent  field  micromechanics  model.  For  the 
viscoelastic  constitutive  equation,  the  cure  and  temperature  dependent  stress  relaxation  modulus  was 
approximated  as  a  Prony  series  of  a  number  of  Maxwell  elements.  The  composite  constitutive  behavior 
is  then  implemented  into  a  general  purpose  finite  element  software  to  determine  the  evolution  of  the 
residual  stress  and  deformations  with  temperature  and  degree-of-cure  dependent  material  behavior. 
Several  numerical  simulations  are  performed  to  compare  the  results  with  other  numerical  and 
approximate  solutions  available  in  the  literature.  Thus,  a  comprehensive  methodology  for  predicting 
process-induced  stress  and  deformations  in  composite  materials  is  formulated,  implemented  and 
validated. 


Why  It  Is  Important: 

Processing  history  strongly  contributes  to  the  residual  stress  development  leading  to  part  warpage  and 
shrinkage  and  potential  reduction  in  the  absolute  strength  of  the  composite.  Several  applications, 
including  modeling  the  residual  stress  and  deformation  evolution  dxiring  RTM  processing,  analysis  and 
optimization  oTprocess  cycles  for  reduced  process  induced  stress,  and  a  novel  use  in  the  prediction  of 
long  term  performance  of  the  composites,  are  anticipated  from  the  developed  methodology. 


Symmetry  Face 


Three  dimensional  geometry  of  a 
typical  part  considered  for  numerical 
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Abstract 


The  process-induced  stress  and  deformation  development  in  unidirectional 
composites  was  investigated.  Processing  history  strongly  contributes  to  the  resid¬ 
ual  stress  development  leading  to  part  warpage  and  shrinkage  necessitating  the 
ability  to  accurately  simulate  the  manufacturing  process  of  composites.  The  fo¬ 
cus  of  this  study  is  the  development  of  a  three-dimensional  thermo-chemical  cure 
simulation  of  the  Resin  Transfer  Molding  (RTM)  process  coupled  with  a  process- 
induced  stress  and  deformation  analysis.  A  Galerkin  finite  element  approach  was 
used  to  solve  the  thermal  equilibrium  and  chemical  kinetics  during  the  curing 
phase  of  RTM  process.  The  part  was  subjected  to  mold  temperature  history 
corresponding  to  a  specified  manufacturing  process  plan.  The  temperature  and 
degree  of  cure  fields  and  their  gradients  axe  obtained  during  the  curing  process. 
The  temperature  and  degree  of  cure  fields  obtained  from  the  cure  simulation  are 
used  to  examine  the  residual  stress  and  deformation  developed  during  the  cur¬ 
ing  phase.  Two  material  constitutive  models,  thermo-elastic  and  viscoelastic, 
were  studied.  For  the  thermo-elastic  approach,  the  rule  of  mixtures  was  used  to 
depict  the  kinetic-viscoelastic  behavior  of  the  resin  modulus.  Effective  compos¬ 
ite  mechanical  properties  were  computed  using  the  instantaneous  resin  and  fiber 
properties  in  a  self-consistent  field  micromechanics  model.  For  the  viscoelastic 
constitutive  equation,  the  cure  and  temperature  dependent  stress  relaxation  mod¬ 
ulus  was  approximated  as  a  Prony  series  of  a  number  of  Maxwell  elements.  The 
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residual  stress  resulting  from  the  processing  history,  thermal  expansion,  and  chem¬ 
ical  shrinkage  were  studied.  For  both  models,  ABAQUS  finite  element  software 
was  implemented  to  perform  the  stress  and  deformation  analysis. 
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NOMENCLATURE 


[A]  coefficient  matrix 

{5}  right-hand-side  vector  of  the  linear  system 

[C\  the  capacity  matrix 

Cpf  specific  heat  of  fiber 

Cpr  specific  heat  of  resin 

Hr  resin  heat  of  reaction 

[K]  the  stiffiiess  matrix 

Kxx  thermal  conductivity  in  x-direction 

Kyy  thermal  conductivity  in  y-direction 

Kgz  thermal  conductivity  in  z-direction 

kf  thefmal  conductivity  of  fiber 

kr  thermal  conductivity  of  resin 

Ne  total  number  of  elements 

Nn  total  number  of  nodes 

[7^]  the  residual  vector  of  linear  system 

R  ideal  gas  law  constant 

it  rate  of  heat  generation 

Ta  reaction  rate 

T  temperature  distribution 

t  time 

{T}  temperature  solution  vector 
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a 

degree  of  cure 

1 

r 

boundary  of  the  mold 

1 

At 

time  step  size  ' 

Q 

time  integration  coefficient 

1 
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Vp 

porosity 

Vf  = 

1  —Vp  fiber  volume  fi-action 

\ 

Pf 

density  of  fiber 

Pr 

density  of  resin 

r 

m 

interpolation  function  vector 

n 

interpolation  function  at  the  ith  node 

n 

the  computational  domain 

1 

the  domaun  of  mold 
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Qp 

the  domain  of  part 
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d/dn 

the  outward  normal  derivative  to  F 

Em 
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instantameous  matrix  modulus 

fully  uncured  matrix  modulus 
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fiilly  uncured  matrix  modulus 
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matrix  modulus  fitting  paurauneter 
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VcH 

volumetric  shrinkage 

*  1' 

totaJ  volumetric  shrinkage  of  matrix 

■  1 
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time-dependent  non-mechamicaJ  straun 
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£eh 

stradn  due  to  chemicail  shrinkage 

Kjk 

effective  chemicail  shrinkage  coefficient 

<l>jk 

effective  thermal  expansion  coefficient 

.  i 

o  time-dependent  stress 
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Q  time-dependent  eflfective  stiffiiess  matrix  of  composite 
Q”  fully  unrelaxed  (elastic)  stiflEness  matrix  of  composite 
Q°°  fully  relaxed  stiffiiess  matrix  of  composite 
Rf  partition  factor 
t,  If  present  and  past  time 

present  and  past  reduced  time 
Or  shift  function 

Tg  glass  transition  temperature 

peak  stress  relaxation  time 
Tjj  discrete  stress  relaxation  time 
Wu  weight  factor 
X  location  vector 
Superscripts 

e  element 

i,  j  the  ith  and  jth  element 
n  the  nth  time  step 


Chapter  1 


INTRODUCTION 


The  use  of  composite  materials  offers  the  attractive  quality  of  high  strength- 
ratios.  This  characteristic  is  beneficial  especially  in  applications  such  as 
the  tail  or  wing  of  an  aircraft.  Although,  composite  materials  offer  tremendous 
advantages,  the  use  of  composite  structures  in  many  applications  has  not  been 
widespread  and  has  been  limited  to  small  components.  This  lack  of  use  can  be 
attributed  to  the  complexify  involved  in  accurately  predicting  the  behavior  of  the 
composite  part.  The  manufacturing  process  has  a  significant  effect  on  the  overall 
mechanical  properties  of  the  part.  In  addition,  the  phases  of  the  manufacturing 
process  are  tightly  coupled, adding  to  the  complexity  of  the  manufacturing  process 
model.  During  the  manufacturing  process,  process-induced  residual  stress  and 
deformation  develop,  which  leads  to  detrimental  effects  such  as  part  warpage  and 
dimensional  instability. 

The  objective  of  this  study  is  to  determine  the  residual  stress  and  defor¬ 
mation  states  for  fiber-reinforced  composite  structures  arising  from  the  manu¬ 
facturing  process,  particularly  the  curing  phase.  During  the  curing  phase,  the 
mechanisms  of  thermal  expansion  and  chemical  shrinkage  drive  the  evolution  of 
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residual  stress  and  defonnation.  In  addition,  the  material  behavior  is  viscoelastic 
necessitating  the  need  of  an  accurate  constitutive  relation  in  terms  of  the  process 
variables  of  time,  temperature,  and  degree  of  cure.  In  order  to  determine  the 
residual  stress,  the  following  approach  was  utilized: 

•  Simulate  the  thermo-kinetic  cure  cycle  process 

•  Obtain  temperature  and  degree  of  cure  fields  during  process  cycle 

•  Model  fiber  reinforcement  structure  to  obtain  effective  mechanical  prop¬ 
erties 

•  Implement- a  time,  temperature,  and  degree  of  cure  dependent  viscoelastic 
constitutive  model  in  a  generalizable  stress  simulation  and  obtain  residual 
stress  and  deformation 

The  methodology  of  the  present  study  is  siunmarized  in  Figure  1.1.  In  this 
methodology,  the  simulation  of  the  cure  phase  and  stress  development  are  not 
performed  simultaneously.'  The  temperature  and  cure  fields  obtained  from  the 
cure  simulation  are  inputs  to  the  residual  stress  simulation. 

The  following  study  will  review  the  modeling  and  simulation  of  residual 
stress  and  deformation  in  polymer  matrix  composites  during  the  cure  cycle.  Chap¬ 
ter  2  focuses  on  the  development  and  implementation  of  the  cure  simulation,  in 
which  the  Galerkin  finite  element  method  is  employed.  NumericaJ  examples  for  a 
flat  plate  and  curve-shaped  geometry  are  presented.  In  Chapter  3,  a  discussion  on 
the  residual  stress  and  deformation  is  presented,  which  includes  a  thermo-elastic 
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and  viscoelastic  formulation.  For  the  viscoelastic  formulation,  a  study  of  the  be¬ 
havior  of  an  isotropic  and  orthotropic  material  was  conducted.  In  addition,  a 
time,  temperature,  and  cure  dependent  viscoelastic  constitutive  relation  was  in¬ 
vestigated.  Numerical  results  of  the  three  different  approaches  Me  also  presented 
after  each  section.  Chapter  4  is  devoted  to  a  discussion  of  the  results. 
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Figure  1.1;  Methodology  of  Coupled  Cure  and  Stress  Simulation 


Chapter  2 


CURE  SIMULATION 


2.1  Introduction 

The  utilization  of  polymer  composite  materials  in  many  of  today’s  indus¬ 
trial,  space,  and  military  applications  has  brought  the  imderstanding  of  the  vari¬ 
ous  liquid  composite  molding  (LCM)  processes  to  the  forefront  of  engineering  and 
science.  In  addition  to  high  strength-to-weight  ratios,  the  use  of  composite  ma¬ 
terials  offer  an  increased  fatigue  life,  improved  corrosion  and  wear  resistance,  low 
thermal  conductivity,  and  design  flexibility.  Unfortunately,  this  design  flexibil¬ 
ity  introduces  greater  design  challenges  not  seen  with  conventional  homogeneous 
materials.  The  various  combinations  of  fiber  and  matrix  configmration,  manu¬ 
facturing  process  and  conditions,  etc.  enhance  the  complexity  of  characterizing 
its  mechanical  properties.  Several  models  have  been  developed  to  determine  the 
mechanical  properties  based  on  the  constituents  and  fiber  volume  fraction.  More¬ 
over,  the  input  parameters  of  the  LCM  process  greatly  affect  the  performance 
and  quality  of  the  finaJ  part.  As  a  result,  the  process  presents  the  manufactur¬ 
ing  challenge  of  producing  a  high  quality  product  at  minimum  cost.  Because  of 
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the  complexity  in  the  chaxacterization  of  composites,  computer  simulations  of  the 
process  are  necessary  to  reduce  the  time  and  cost  of  experiments  and  trial  and 
error  manufacturing. 

One  of  the  typical  LCM  processes  used  to  develop  composites  is  resin  trans¬ 
fer  molding  (RTM).  RTM  is  often  used  to  produce  large,  net-shaped  parts  yielding 
high  stfen^h-to-weight  ratios.  The  process  entails  a  thermosetting  resin  injected 
into  a  heated  mold  packed  with  a  fiber  preform  and  cured  under  controlled  tem¬ 
peratures.  The  process  cycle  for  polymer  matrix  composites  (PMCs)  with  a  ther¬ 
mosetting  resin  matrix  consists  of  the  following  stages:  1)  mold  filling,  2)  curing, 
3)  mold  cooling,  .^d  4)  demolding.  During  the  mold  filling  stage,  the  injection 
pressure,  temperature  of  the  mold,  permeability  of  the  matrix,  resin  viscosity,  flow 
rate,  location  of  gates,  etc.  dictate  the  quality  of  the  part.  While  this  stage  of  the 
process  plays  .an  important  role  in  the  outcome  of  the  final  part,  the  mold  filling 
and  curing  stages  are  decoupled  and  the  curing  stage  is  the  main  focus  of  this 
study. 

A  typical  cure  cycle  consists  of  two  dwell  temperatures.  In  a  two-step  cure 
cycle,  the  material  temperature  is  raised  from  room  temperature  to  the  first  dwell 
temperature  and  held  constant  for  about  an  hour.  During  the  first  dwell,  excess 
gases  are  removed  and  the  viscosity  of  the  matrix  material  is  lowered  to  facilitate 
flow  and  compaction  of  the  part.  The  temperature  is  increased  to  the  second 
dwell  temperature  which  allows  the  cross-linking  of  the  poljnner  to  take  place. 
[1]  There  is  a  minimum  temperature  required  before  cross-linking  begins.  High 
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cure  temperatures  increase  the  reaction  rate  and  lower  the  cycle  time  and  cost. 
However,  higher  cure  temperatures  lead  to  the  development  of  process-induced 
residual  stresses,  which  result  in  part  warpage  and  dimensional  instability.  There¬ 
fore,  it  is  important  to  understand  and  simulate  the  mechanisms  contributing  to 
the  curing  phase  stress  development. 

-A distinction  between  thin  and  thick  section  composites  must  be  made  be¬ 
fore  proceeding  with  any  type  of  residual  stress  analysis.  Thick  section  composites 
experience  a  number  of  manufacturing  problems  that  do  not  appear  in  thin  section 
laminates  such  as  thermal  spiking,  non-uniform  consolidation,  and  non-uniform 
curing  [2].  The  heat  released  from  the  exothermic  cross-linking  reaction  causes 
the  temperature  increase  in  the  interior  of  the  part.  With  thin  section  composites, 
the  heat  is  quickly  transmitted  to  the  outer  edges.  However,  in  thick  section  com¬ 
posites,  the  low  thermal  conductivity  prevents  the  quick  dissipation  of  the  heat 
resulting  in  thermal  spiking  and  large  thermal  gradients. 

2.2  Literature  Review 

^  Large  thermal  gradients  create  differential  degree  of  cure  and  cause  varying 
thermal  expansions  from  the  mismatch  of  coeflBcients  of  thermal  expansion  (GTE) 
and  chemical  shrinkage,  especially  in  thick  section  composites  [3],  which  lead  to 
significant  residual  stress  development.  Various  studies  have  been  conducted  in 
order  to  gain  an  understanding  of  the  curing  phase.  Early  work  by  Loos  and 
Springer  [4]  involved  the  development  of  a  one-dimensional  cure  simulation  of  a 
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flat  plate  using  an  implicit  finite  difference  method.  Bogetti  and  Gillespie  [5] 
conducted  a  two-dimensional  cure  simulation  in  thick  anisotropic  thermosetting 
composites  with  a  coupled  thermal  and  chemical  kinetics  formulation  using  the 
finite  difference  method.  Yi  et.  al.  [6]  developed  a  nonlinear  transient  heat 
transfer  finite  element  model  where  material  parameters  are  both  temperature 
and  de^ee  of  cure  dependent.  In  the  work  of  Park  and  Lee  [7],  a  two  dimensional 
cure  simulation  was  also  developed  and  introduced  a  degeneration  method  to  build 
the  thermal  conductivity  matrix.  Zhu  et.  al  [8]  presented  a  full  three-dimensional 
coupled  thermo-chemo-viscoelastic  model  simulating  the  heat  transfer,  curing, 
residual  stresses  and  deformation  of  a  composite  part.  For  the  cure  simulation, 
the  thermal  and  cure  mass  matrices  were  Imnped  via  the  HRZ  lumping  scheme 
[9]  to  prevent  negative  degree  of  cure  values.  The  Newton-Raphson  method  and 
an  adaptive  time  stepping  technique  are  employed  in  the  finite  element  solution. 

The  thermo-chemo  equilibriiim  equation  includes  the  internal  heat  genera¬ 
tion  term  resulting  from  the  exothermic  cure  reaction.  Diflferent  reaction  rate  mod¬ 
els  have  been  presented  for  several  material  systems,  particularly  glass/polyester 
and  graphite/epoxy.  The  model  for  the  Hercules  3501-6  epoxy  resin  was  developed 
by  Loos  and  Springer  [4].  Bogetti  zind  Gillespie  [5]  modeled  the  glass/polyester  sys¬ 
tem.  Later  cure  models  [l],[6]-[8],  and  [10]  used  the  same  models  for  the  polyester 
and  epoxy  systems. 
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2.3  Problem  Description 

In  a  study  of  the  residual  stress  and  deformations  of  PMCs,  the  temperature 
and  degree  of  cure  history  throughout  the  part  is  necessary.  The  Galerkin  finite 
element  formulation  was  used  for  the  coupled  thermo-kinetic  simulation  in  this 
effort.  This  methodology  is  suited  for  solving  large  meshes  using  element-free  and 
iterative  solution  tediniques  for  increased  computational  efficiency  when  dealing 
with  large  mesh  sizes  [11].  The  paper  presents  the  formulation,  illustration  of 
typical  results,  and  examples  used  for  verification  of  the  implementation.  Two 
numerical  examples  are  presented  for  which  prior  results  and  experimental  data 
are  available  in  the'  literature.  The  next  section  presents  the  thermo-kinetic  for¬ 
mulation  proceeded  by  the  two  numerical  examples.  The  cure  simulation  is  then 
followed  by  different  approaches  to  model  the  residual  stress^  and  deformations 
during  the  curing  phase 

2.4  Cure  Simulation  Formulation 

The  in-mold  curing  phase  of  the  RTM  process  is  formulated  as  a  coupled 
thermo-kinetic  analysis.  The  presence  of  separate  and  dissimilar  material  phases 
(fiber  and  matrix)  is  modeled  using  effective  property  homogenization.  The  ef¬ 
fective  physical  and  thermal  properties  are  determined  from  the  fiber  and  matrix 
properties,  the  fiber  volume  fraction  and  orientation.  Two  separate  kinetics  mod¬ 
els  are  used  to  simulate  the  cure  reaction,  one  appropriate  for  the  glass/polyester 
and  the  other  for  carbon/epoxy  systems.  The  temperature  cycle  used  for  the 
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curing  of  part  is  tramslated  into  time-dependent  mold  boundary  conditions. 

2.4.1  Energy  Balance 

The  three-dimensional  time  dependent  energy  balance  with  orthotropic  ma¬ 
terial  conductivity  is  given  as  follows: 

-  .  Qjn  q2J’  ^rp 

{VpPrCpr  +  Qy2  dz^  ^  (^-l) 

where  R  is  the  rate  of  heat  generation  from  the  chemical  reaction,  which  is  related 
to  the  total  heat  of  reaction  {Hr)  and  rate  of  reaction  (ra)  as  given  by  the  following 
equation: 

R  =  Hrra  (2.2) 

The  equivalent  thermal  conductivities,  Kxx,  Kyy  and  iTzz,  are  computed  from  the 
fiber  and  matrix  properties.  Kxx  is  calculated  from  rule  of  mixtures  [12]  using  the 
fiber  volume  fraction  Vf.  The  transverse  (Kyy)  eind  through- the-thickness  (Kzz) 
conductivity  values  au:e  similaurly  calculated  based  on  the  fiber  orientation  details. 


^  kr  —  kf 

kr  -{■  kf  Vf{kr  kf) 

(2.3) 

=  Vfkf{l  +  Vpkt)  -1-  Vpkr{l  -  Vfkt) 

(2.4) 

The  governing  equation  is  solved  subjected  to  the  mold-wall  temperature  bound¬ 
ary  conditions  and  vanishing  heat  flux  conditions  on  symmetry  boundaries; 

T  =  on  mold  surfaces  (2.5) 

—K—  =  0  on  symmetry  surfaces 
on 
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2.4.2  Cure  Kinetics 

Curing  of  polymers  is  temperature  dependent  as  the  reaction  kinetics  vary 
with  temperature.  Also,  the  reactions  are  exothermic  hence  the  coupling  between 
thermal  balance  and  cure  kinetics.  The  reaction  rate  and  the  degree  of  cure,  a 
are  related  by  the  following  rate  equation: 

da 

a  =  ’■"  (2-6) 

with  the  initial  condition,  Oq  =  ai  at  t  =  U. 

For  different  types  of  polymerization  appropriate  kinetic  models  are  required 
to  determine  the. rate  and  heat  of  reaction.  Two  different  models  were  used  for 
modeling  the  glass/epoxy  and  graphite/epoxy  systems.  The  kinetic  model  used 
for  the  glass/polyester  material  is  as  follows  [5]: 


^  =  (A-,+ir2a”)(l-a)" 

(2.7) 

Ki  =  Alt 

(2.8) 

,  K2  = 

The  model  for  the  graphite/ epoxy  is  a  three-term  rate  equation  as  shown  below 
[4]: 


dct 

—  =  {Ki  +  A'2a)(l  -  Q!)(0.47  -  a)  for  (a  <  0.3) 

(2.9) 

da  ^ 

, 

=  ^3(1  -  tt)  for  (a  >  0.3) 

(2.10) 

Ki  =  Aie  kt 

(2.11) 

K2  =  A2e  RT 
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Ks  =  Ase  RT 

The  thermal  balance  and  kinetics  equations  are  solved  in  an  incremental/iterative 
time  stepping  scheme  where  the  temperature  and  degree  of  cure  is  determined 
sequentially  at  each  time  step. 


2.4:3 '  Galerkin  Finite  Element  Method 


The  computational  domain,  is  discretized  by  a  union  of  finite  elements 
as  Q  =  S  U^i  where  D  =  0,  if  i  ^  j  and  Nf,  is  the  total  number 
of  elements  in  f2.  The  approximate  element  solution  is  defined  as: 

T=^{^x,y,z)fmt)}  (2.12) 

For  the  Galerkin  method,  the  weight  function  is  chosen  as: 


w  =  {$(a:,  y,  z)} 


(2.13) 


Equation  (2.1)  can  be  integrated  by  multiplying  the  weight  function  and  using 
Green’s  theorem, 

a/'ri 

(2.14) 


Id  ^ + w  m = m 


dt 


where 


ld  =  EX„.WW’’''*^ 


(2.15) 


e 

\T 


dx  dx  dy  dy  dz  dz 


Ck  = 


’Ti 


VpPr^pr  "b  "^fPf^pf 


dfl 
(2.16) 
(2.17) 


I 
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(2.18) 

1 

-r  ^ 

Q  _  Pr'^pHr 

VpPrCpr  +  VfpfCpf 

(2.19) 

K 

1 

Equation  (2.14)  is  discretized  by  the  backward  Euler  method  in  time  where  At  = 

”  1 

{|C)  +  At  M)  {■n'"+‘>  =  10  +  At  {R} 

(2.20) 

Equation  (2.20)  can  be  simply  denoted  as  a  system  of  linear  algebraic  equations 

defined  by, 

'1 

J 

(2.21) 

^  1 

where, 

f 

[>4]  =  [C]  +  At[X:] 

(2.22) 

{B}  =  [c]{r}^^^  +  At{n} 

(2.23) 

By  multiplying  by  the  weight  function,  Eq.  (2.6)  can 

be  integrated  in 

the  same  fashion  as  the  thermal  balance  equation. 

■  =  w 

(2.24) 

1 

The  kinetics  equation(2.24)  is  also  discretized,  where  At  =  t("+^) 

- 1(”). 

[C] =  [C]a(”)  + At  TO 

(2.25) 

where 

-i 

(2.26) 

•  -« 

.  i 

(2.27) 
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Equation  (2.25)  can  be  denoted  as  a  system  of  linear  algebraic  equations  defined 

by, 

=  (2.28) 

where, 

{B}  =  [C]  +  At  {7^}  (2.29) 

HRZ  mass  lumping  scheme  [8]  is  employed  to  diagonalize  the  cure  matrix,  [C\. 
The  diagonalization  leads  to  a  inexpensive  solution  of  Eq.  2.29  at  each  increment. 

The  solution  is  performed  with  a  banded  solver  using  Gibbs-Poole-Stockmeyer 
bandwidth  (ACM  CALGO  Algorithm:582)  reduction.  The  computational  times 
recorded  for  a  363  degree  of  freedom  problem  are  154  milliseconds  per  time  step  on 
a  500  MHZ  Pentium  III  class  proc^sor.  For  a  problem  with  a  larger  mesh  10125 
degrees  of  freedom  and  initial  band  width  of  5343,  the  bandwidth  reduction  al¬ 
gorithm  reduced  the  bandwidth  to  287,  and  the  computational  time  required  was 
700  milliseconds  per  .time  step. 

2.5  Numerical  Examples 

Two  numerical  examples  are  presented  to  illustrate  the  results  obtained. 
The  results  from  the  current  method  are  compared  with  those  avedlable  in  the  liter¬ 
ature  and  experimental  data  when  available.  Two  geometry  choices,  namely  a  flat 
plate  and  a  curve-shaped  part,  and  two  material  systems,  namely  glass/polyester 
and  carbon/epoxy,  are  used  for  the  numerical  examples. 
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Table  2.1:  Material  Properties  of  Glass/Polyester 


Density  (p) 

1890  kg/m* 

Polymer  specific  heat  (Cp) 

1260  J/kg-K 

Thermal  Conductivity 

0.2163  W/m-K 

Kxx/^ZZ 

2 

2.5.1  Curing  of  a  Flat  Plate 

-  A  flat  plate  with  dimensions  of  15.24  cm  x  15.24  cm  and  a  uniform  thickness 
of  2.54  cm  is  considered.  Due  to  symmetry,  a  quarter  of  the  plate  is  meshed.  The 
geometry  of  the  model  analyzed  was  7.62  cm  x  7.62  cm  and  1.27  cm  thick.  The 
thermal  and  kinetic  properties  for  the  glass/polyester  can  be  found  in  Tables  2.1 
and  2.2,  respectively.  The  flnite  element  model  used  composed  of  a  200  element 
8  node  brick  mesh.  The  part  was  subjected  to  a  prescribed  temperature  cycle  for 
which  experunental  results  were  available  [5].  The  cycle  entails  a  60  minute  dwell 
at  78®C,  a  rise  to  OO^C  and  a  70  minute  dwell  at  126‘’C  as  seen  in  Figure  2.1. 
Figures  2.2  and  2.3  show  the  temperature  and  cure  evolution  at  the  center  of  the 
plate  for  the  glass/polyester  system.  These  results  are  in  good  agreement  with 
the  published  experimental  results.  Figure  2.2  shows  a  maximum  temperature  of 
128®C  (a  thermal  spike)  occurring  at  165  minutes  into  the  cure  cycle  due  to  the 
exothermal  curing  reaction.  There  was  a  38°C  temperature  difference  between  the 
boundary  and  the  center  of  the  part.  The  accuracy  of  the  maximum  temperature 
was  found  to  be  sensitive  to  the  time  step  utilized  in  the  numerical  analysis.  A 
time  increment  of  1  sec  was  adequate  to  capture  the  effects  of  the  exotherm  in 
the  simulation.  The  simulation  of  the  glass/polyester  plate  was  completed  with 
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Table  2.2:  Cure  Kinetics  for  Glass/Polyester 


Ai  (min~l) 
A2  (min~l) 
Bi  (J/mole) 
E2  (J/mole) 
m 
n 

R  (J/mol-K) 
ffr  (J/kg) 


0 

3.7  X  10^2 
0 

1.674  X  10^ 
0.524 
1.476 
8.31434 
77,500 


13,920  increments  (4  hours  of  simulation  time)  with  a  CPU  time  of  2146  seconds 
(36  minutes)  on  a  personal  computer.  Prom  the  plot  of  the  degree  of  cure  in  Figure 
2.3,  the  value  of  the  degree  of  cure  goes  from  0  to  0.99  in  a  short  period  of  time 
from  130-150  minutes.  This  rapid  change  occurs  during  the  period  of  the  thermal 
spike  further  emphasizing  the  sensitivity  of  the  simulation  to  the  time  step. 


0  50  100  '  150  200  250  300  350  400 

TIME  (MIN) 


Figure  2.1:  Boundary  Temperature  Cycle  for  Glass/Polyester  Material 


21 


0  50  100  150  200 

-  TIME  (MIN) 

Figure  2.2:  Temperature  Profile  for  Center-point  Node  for  Glass/Polyester  Plate 
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The  plate  curing  simulation  was  also  run  with  graphite/epoxy  material  prop¬ 
erties  as  given  in  Tables  2.3  and  2.4.  The  temperature  cycle  for  the  graphite/epoxy 
material  [13]  found  in  Figure  2.4  had  a  1  hour  dwell  at  116®C  and  a  2  hour  dwell 
at  177"C.  The  graphite/epoxy  plate  results  in  Figures  2.5  and  2.6  were  also  in 
agreement  with  the  experiment  in  [13].  Because  of  the  change  in  cure  kinetic 
models'when  the  degree  of  cure  reaches  above  0.3,  the  degree  of  cure  profile  ap¬ 
pears  to  have  two  peaks  with  the  transition  at  104  minutes.  There  are  two  thermal 
spikes  in  the  temperature  profile  of  the  center  point  node  in  Figure  2.5.  For  the 
first  spike,  there  was  a  O^C  difference  compared  to  the  boundary.  There  was  a 
31®C  difference  for  the  second  spike.  The  cure  simulation  for  the  graphite/epoxy 
plate  required  10,000  iterations  (26  minutes  of  CPU  time)  with  a  1  sec  time  step. 
Although  the  flat  plate  geometry  with  uniform  thickness  may  not  warrant  a  three- 
dimensional  analysis,  this  case  was  used  to  benchmark  and  verify  the  performance 
of  the  present  method. 
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Table  2,3:  Material  Properties  of  Graphite/Epoxy 


Density  (p) 

1578  kg/m^ 

Polymer  specific  heat  {Cp) 

862  J/kg-K 

Thermal  Conductivity 

0.4135  W/m-K 

J^xxf  J^ZZ 

1,  5, 10 

Table  2.4:  Cure  Kinetics  for  Graphite/Epoxy 


Ai  {min~l) 

2.102  X  10» 

'A2  {mirrl) 

-2.014  X  10® 

A3  {min~l) 

1.960  X  10® 

El  (J/mole) 

8.07  X  10^ 

E2  (J/mole) 

7.78  X  10^ 

E3  (J/mole) 

5.66  X  10^ 

R  (J/mol-K) 

8.31434 

Hr  (J/kg) 

198.6  X  10® 

TEMPERATURE  (DEQC) 


Figure  2.5:  Temperature  Profile  for  Center-point  Node  for  Graphite/Epoxy  Plate 
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2.5.2  Curve-Shaped  Part 

As  the  second  example,  the  curve-shaped  part  presented  in  [8]  is  analyzed. 
The  dimensions  of  the  cross  section  are  shown  in  Figure  2.7.  The  part  had  an 
overall  width  of  8  cm.  Because  of  half-symmetry,  the  model  had  a  width  of  4  cm. 
The  material  system  used  for  the  part  is  graphite/epoxy  with  properties  shown 
in  Tables  3  and  4.  The  temperature  cycle  in  Figure  2.4  was  applied  to  all  the 
boundaries  except  the  front  face  because  of  symmetry.  The  mesh  employed  has 
448  elements  and  725  nodes  with  a  bandwidth  of  34.  The  temperature  and  degree 
of  cure  profiles  for  a  node  in  the  center  of  the  part  can  be  fotmd  in  Figures  2.8  and 
2.9,  respectively.  Using  a  1  sec  time  step,  the  simulation  was  complete  after  9,960 
iterations.  The  computational  time  required  was  in  the  order  of  130  milliseconds 
per  iteration. 

For  this  example  we  present  a  series  of  results  detailing  temporal  and  spatial 
variation  of  the  temperature  and  degree  of  cure  within  the  part.  Three  points  in 
time  are  selected  for  obtaining  the  spatial  variations  of  temperature  and  cure  in 
the  part.  The  first  point,  ta,  is  at  100  minutes  into  the  cure  cycle,  before  the 
thermal  rise.  The  other  two  are  at  ti,=  113  min  and  tc  -  125  minutes.  The 
final  time,  tc,  is  at  the  peak  of  the  thermal  spike,  while  the  second  point  is  the 
midpoint.  These  points  correspond  to  the  rapid  reaction  phase  of  the  cure  cycle 
as  seen  in  Figure  2.8.  The  temperature  and  degree  of  cure  contours  are  shown  on 
four  separate  cut-planes  in  the  part  seen  in  Figures  2.10-2.13.  As  seen  in  these 
figures,  the  temperature  increases  from  a  low  temperature  of  116®C  to  a  high 
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temperature  of  195"C  within  a  short  period  of  time,  approximately  25  minutes. 
The  temperature  difference  is  79°C.  Also  in  this  time  period,  the  degree  of  cure 
rises  from  0.276  to  0.699,  a  difference  of  0.41.  The  contour  plots  also  show  that 
the  highest  temperature  and  degree  of  cure  occur  mainly  in  the  middle  portion  of 
the  part.  At  point  1  and  2,  the  temperature  and  degree  of  cure  is  mainly  uniform 
throughout  the  part,  whereas  point‘3  displays  a  temperature  difference  of  18"C. 


TIME  (MIN) 

Figure  2.8:  Temperature  Profile  for  Center-point  Node  for  Curve  Shaped  Part 
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Figure  2.10:  Curved  Section:  Temperature  (left)  and  Degree  of  Cure  (right)  Con¬ 
tours 


time  =  100  min 


34 


CVREJOOO 

A2M 


time  =  113  min 


i’"  ^ 


(a) 

(d) 

1 

J 

ctmsjm 


time  =  125  min 


CUMJSOO 

AMf 


Figure  2.11.  Vertical  Section:  Temperature  (left)  and  Degree  of  Cure  (right) 
Contours 


Chapter  3 


RESIDUAL  STRESS  AND  DEFORMATION 


3.1  Introduction 

Composite  materials  offer  a  number  of  desirable  qualities  as  previously 
stated.  However;  their  use  has  been  limited  by  the  complex  issues  involved  in 
manufacturing  a  dimensionally  accurate  part.  The  residual  stress  and  deforma¬ 
tion  development  during  the  curing  phase  has  a  significant  effect  on  both  the 
dimensions  and  mechanical  performance  of  the  part.  Residual  stresses  can  result 
in  warpage,  matrix  cracks,  and  delamination  causing  part  failure  earlier  than  ex¬ 
pected.  The  ability  to  predict  the  behavior  of  composite  materials  is  necessary 
especially  when  being  used  in  applications  such  as  aircrafts  and  space  shuttles. 
The  predictability  of  the  composite  strength  would  aid  in  the  the  design  of  com¬ 
posite  parts  with  better  dimensional  stability.  In  addition,  an  accurate  prediction 
of  composite  behavior  would  eliminate  the  expenses  involved  in  manufacturing 
by  trial  and  error.  Therefore,  understanding  the  mechanisms  involved  in  residual 
stress  development  and  creating  computer  simulations  of  the  process  is  important 
and  this  is  the  focus  of  this  research.  Residual  stresses  mainly  arise  from  thermal 
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expansion  coefficient  mismatch  between  the  composite  constituents,  volumetric 
chemical  shrinkage  of  the  resin,  and  non-uniform  curing.  These  factors  are  espe¬ 
cially  significant  in  the  case  of  thick  section  composites  where  large  thermal  and 
degree  of  cure  gradients  exist.  These  gradients  induce  spatially  varying  thermal 
expansions  and  chemical  hardening  effects  where  there  is  nonuniform  mechanical 
property  development  across  the  laminate  thickness  [13]. 

3.2  Literature  Review 

Several  studies  have  been  conducted  in  the  analysis  of  process  induced  resid¬ 
ual  stresses.  The  earliest  study  of  residual  stresses  in  thermoset  matrix  composites 
was  performed  by  Hahn  and  Pagano  [14].  In  their  study,  a  classical  lamination 
theory  (CLT)  was  applied  to  their  elastic  analysis  of  a  thin  laminate,  in  which 
a  stress-free  state  before  cool-down  was  assumed.  Other  studies  employing  the 
same  technique  simulated  the  in-plane  stresses  and  curvatures  [15]-[17].  Radford 
[18]  also  used  the  classical  laminated  plate  theory  with  volume  fraction  varia¬ 
tions  to  study  the  cure  shrinkage  induced  warpage  in  flat  uni-axial  carbon/epoxy 
composites.  Results  showed  that  the  volume  fraction  gradients  induced  from  the 
manufacturing  process  were  an  important  component  to  part  warpage. 

Early  work  by  Weitsman  [19]  showed  that  the  viscoelastic  response  of  the 
resin  must  be  taken  into  account.  The  residual  stresses  resulting  from  cool-down 
was  reduced  by  twenty  percent  in  comparison  to  a  linear  elastic  analysis.  An  op¬ 
timal  cool  down  path  was  derived  based  on  this  analysis  in  a  follow  up  study[20]. 
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In  the  work  of  Hodges  et  al.  [21],  the  effect  of  cure  temperature  on  the  resid¬ 
ual  stresses  for  an  epoxy  resin  (Cyanamid-Fothergill)was  studied.  A  lower  cure 
temperature  was  proposed  to  reduce  the  mismatch  in  the  thermal  expansions.  A 
number  of  curing  cycles  were  considered  and  monitored  with  measurements  of  vis¬ 
cosity  and  volume  changes  using  differential  scanning  calorimetry  (DSC).  Stango 
and  Wang  [22]  studied  the  viscoelastic  response  of  a  laminate  by  examining  the 
interlaminar  stresses  at  the  free  edge  due  to  cool-down. 

Kim  and  Hahn  [23]  monitored  the  warpage  of  unsymmetric  cross-ply  lami¬ 
nate  and  used  an  elastic  method  to  relate  the  warpage  to  reiduaJ  stress.  Bogetti 
and  Gillespie  [24],  investigated  the  development  of  residual  streses  during  cure 
in  thick  composites  for  a  glass/polyester  material  system.  Their  model  took  the 
matrix  shrinkage  and  thermal  expansion  into  account  by  contributing  to  the  total 
strains.  From  the  incremental  total  strains,  the  stress  and  deformations  were  cal¬ 
culated  based  on  a  cure-dependent  resin  modulus.  Their  work  was  based  on  an 
elastic  model,  which  coupled  a  one-dimensional  cure  simulation  with  laminated 
plate  theory.  This  method  of  determining  residual  stresses  was  later  adopted  in  the 
work  of  Teplinksy  and  Gutman  [10],  in  which  the  effects  of  stacking  geometry  and 
laminate  thickness  on  stresses  and  strains  were  examined.  A  linear,  viscoelastic 
stress  model  was  derived  by  White  and  Hahn  [1]  to  include  the  effects  of  chemi¬ 
cal  and  thermal  strains  during  autoclave  or  hot  press  processing  of  PMCs.  The 
mechanical  properties  had  a  functional  dependence  on  the  degree  of  cure  and  the 
transverse  compliance  was  assumed  to  be  the  only  time  dependent  compliance. 
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The  residual  moments  and  curvatures  for  unsymmetric  cross-ply  laminates  were 
calculated  and  validated  for  a  graphite/bismaleimide  composite  system  in  a  com¬ 
panion  paper  [25].  The  study  showed  that  the  contribution  of  the  chemical  strains 
to  residual  stress  was  less  than  4%.  Unfortunately,  this  model  was  restricted  to 
thin  laminates. 

In  the  work  of  Kim  and  White  [26]  and  [13],  residual  stresses  in  thick 
graphite/epoxy  laminates  were  studied  using  a  viscoelastic  analysis,  where  the 
time-dependent  stiffnesses  were  modeled  using  a  series  of  exponential  functions 
and  based  on  cure  dependent  relaxation  times.  It  was  found  that  the  residual 
stress  developed  in  thick  section  composites  was  significant  especially  in  the  trans¬ 
verse  direction.  The  model  included  the  effects  of  strain  from  thermal,  moisture 
and  chemical  expansions.  The  effect  of  chemical  shrinkage  was  found  to  be  neg¬ 
ligible  for  the  material  system  analyzed.  Also,  an  elastic  model  was  found  to  be 
adequate  if  the  stress  development  is  confined  to  the  cool-down  stage  and  chem¬ 
ical  shrinkage  and  hardening  effects  aren’t  significant.  A  recent  study  using  the 
same  methodology  as  Kim'  and  White  was  conducted  by  Zhu  et  al.  [8].  A  three- 
^mensional  coupled  thermo-chemo-viscoelastic  model  was  developed  to  simulate 
the  heat  transfer,  curing,  residual  stress  and  deformation  during  the  entire  cure 
cycle.  Numerical  results  yielded  significant  stress  development  prior  to  cool-down. 

In  order  to  describe  the  viscoelastic  response  of  composites,  Kim  and  White 
[27]  conducted  a  study  into  the  stress  relaxation  of  3501-6  epoxy  resin  at  cure 
states  ranging  from  0.57  to  0.98.  A  Prony  series  was  used  to  describe  the  relaxation 
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modulus  and  shift  functions  used  to  obtain  reduced  times  were  derived  based  on 
curve  fitted  data.  The  experimental  study  showed  that  the  cure  state  had  a 
profound  effect  on  the  stress  relaxation  of  epoxy.  Further  study  of  the  stress 
relaxation  of  3501-6  epoxy  resin  was  done  by  White  and  Hartman  [28]  in  which 
the  specimens  were  tested  in  three-point  bending  to  obtain  creep  compliance  over 
a  wide  temperature  range.  Direct  inversion  and  the  Hopkins-Hamming  method 
were  used  to  calculate  the  stress  relaxation  modulus  mater  curves  from  the  creep 
compUance  curves.  Chemical  aging  effects  were  also  taken  into  account.  Another 
experimental  method  for  viscoelastic  characterization  was  presented  by  O’Brien 
et  al  [29]  using  parallel  plate  rheometry  to  measure  the  material  behavior  prior  to 
the  gel  point. 

The  development  of  the  mechanical  properties  during  cure  was  also  exam¬ 
ined  by  Simon  et  al.  [30]  by  using  a  methodology  combining  the  principles  of 
time-temperature  superposition,  time-cross-link  density  superposition,  and  the 
rubbery  elasticity  of  networks.  The  shift  factors  in  this  formulation  followed  the 
Williams-Landel-Ferry  (WLF)  equation.  The  viscoelastic  model  was  extended  to 
predict  the  isotropic  residual  stresses  in  a  commercial  epoxy  in  [31].  Like  previous 
viscoelastic  models,  the  shear  modulus  was  modeled  as  a  sum  of  Maxwell  elements 
each  having  a  characteristic  relaxation  time.  The  results  of  the  study  showed  that 
cure  shrinkage  contributes  significantly  to  the  total  residual  stresses  and  an  elastic 
analysis  overestimates  the  residual  stresses.  Chemical  shrinkage  was  found  to  be 
a  major  contributor  to  the  spring-forward  effect  in  L-shaped  continuous  fiber  lam- 
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mates  in  Wiersma  et  al.  [32].  The  prediction  of  the  spring-forward  was  sensitive 

'I 

to  the  amount  of  cure  shrinkage  and  viscosity  profile  during  the  cure  cycle  using  I 

a  viscoelastic  model.  > 

Extensive  studies  in  the  viscoelastic  stress  analysis  of  composites  was  also 
conducted  by  Yi  et  al.  In  [33]  and  [34],  a  finite  element  formulation  for  analyz¬ 
ing  interlaminar  stress  fields  in  nonlinear  amsotropic  viscoelastic  composites  was 
presented.  Schapery’s  single  integral  formulation  was  extended  to  account  for 

'I 

viscoelastic  anisotropy  and  multiaxial  stress  states.  The  finite  element  procedure  ^ 

was  also  used  to  study  the  behavior  during  cool-down  in  [35].  Hilton  and  Yi  [36]  i 

performed  a  stoch^tic  delamination  simulation  of  nonlinear  viscoelastic  compos¬ 
ites.  A  similar  methodology  using  Schapery’s  nonlinear  viscoelastic  model  was 
presented  in  [37]  and  [38]  for  unidirectional  carbon/epoxy  composites..  A  nonlin¬ 
ear  three-diniensional  viscoelastic  analysis  was  performed  for  unidirectional  and 
cross-ply  laminates  by  Chen  et  al.  [39] with  a  viscoelastic  model  represented  by  a 
finite  series  of  Kelvin  elements  coupled  with  an  elastic  spring.  A  representative  | 

volume  element  was  analyzed  rather  than  the  whole  composite  geometry.  The  ; 

showed  that  the  viscous  behavior  of  the  matrix  plays  an  important 
role  in  the  evaluation  of  the  residual  stresses.  Fisher  and  Brinson  [40]  also  per-  ; 

formed  a  micromechanicaJ  investigation  using  the  Mori-Tanaka  (MT)  method  and 
the  Benveniste  solution  to  include  the  effects  of  viscoelastic  interphase  regions. 
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3.3  Problem  Description 

The  contribution  of  the  spatial  and  temporal  variance  of  the  composite  ma¬ 
terial  behavior  during  the  curing  phase  to  the  evolution  of  stress  and  deformation 
is  the  focus  of  the  present  study.  Most  of  the  models  on  process-induced  stress 
from  previous  studies  are  limited  to  the  cool-down  phase  of  the  composite.  In  ad¬ 
dition,  several  finite  element  codes  were  developed.  However,  presently,  a  generic 
platform  to  simulate  the  manufacturing  process  does  not  exist.  The  objective  is 
to  study  the  behavior  of  the  composite  and  evolution  of  residual  stress  and  defor¬ 
mation  during  the  curing  phase  to  aid  in  the  development  of  a  generic  architecture 
to  simulate  the  coupled  phases  of  the  composite  manufacturing  process. 

The  material  constitutive  models  adopted  are  coupled  with  the  three-dimensional 
cure  simulation  to  imderstand  the  mechanisms  involved  in  the  stress  develoi>- 
ment.  The  models  take  into  account  the  effects  of  thermal  expansion,  chemical 
cure  shrinkage,  and  non-uniform  curing,  which  are  all  important  sources  of  in¬ 
ternal  loading.  The  graphjte/epoxy  composite  system  was  used  to  illustrate  the 
significant  stress  development  during  the  curing  phase;  thereby  emphasizing  the 
need  for  a  stress  analysis  prior  to  cool-down.  Two  material  constitutive  models, 
thermo-elastic  and  viscoelastic,  were  studied.  For  the  thermo-elastic  approach, 
the  rule  of  mixtures  was  used  to  depict  the  kinetic-viscoelastic  behavior  of  the 
resin  modulus.  Effective  composite  mechanical  properties  were  computed  using 
the  instantaneous  resin  and  fiber  properties  in  a  self-consistent  field  micromechan¬ 
ics  model.  For  the  viscoelastic  constitutive  equation,  the  cure  and  temperature 
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dependent  stress  relaxation  modulus  was  approximated  as  a  Prony  series  of  a 
number  of  Maxwell  elements.  For  both  models,  ABAQUS  finite  element  software 
was  implemented  to  perform  the  stress  and  deformation  analysis  on  a  composite 


system  of  unidirectional  fibers. 
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3.4  Thermo-Elastic  Formulation 

During  the  curing  process,  the  state  of  the  thermosetting  resin  dictates  the 
effective  mechanical  properties  of  the  composite  system.  The  modulus  of  the  resin 
goes  from  low  to  high  as  the  resin  changes  from  a  liquid  viscous  state  to  a  vis¬ 
coelastic/elastic  solid.  Prom  Figure  3.1,  the  evolution  of  the  resin  modulus  can  be 
separated  into  three  regions  [41].  At  the  beginning  of  the  cure  cycle,  the  resin  is 
assumed  to  be  a  fully  uncured,  viscous  fluid  where  the  stiffiiess  is  negligible.  In  the 
second  region,  significant  increases  in  the  stiffness  and  chemical  shrinkage  occur 
due  to  the  cross-linking  reaction.  During  this  phase,  the  resin  quickly  goes  from 
an  uncured  to  a  cured  state.  Also  in  this  region,  the  chemical  kinetic  hardening 
and  viscoelastic  relaxation  are  competing  mechanisms  governing  the  mechanical 
properties  of  the  resin.  The  behavior  of  the  resin  in  the  third  region  is  observed 
to  be  viscoelastic  at  high  temperatures  and  elastic  at  lower  temperatures.  For 
the  thermo-elastic  formulation,  the  properties  of  the  fiber  were  assumed  to  be 
constant  and  independent  of  cure.  The  instantaneous  resin  modulus  was  modeled 
as  a  function  of  degree  of  cure  and  used  in  conjunction  with  the  constant  fiber 
properties  and  volume  fraction  in  a  micromechanics  model  to  determine  the  ef¬ 
fective  composite  properties.  While  the  resin  modulus  was  strongly  influenced  by 
the  curing  process,  the  resin  Poisson  ratio,  was  assumed  to  be  constant.  In 
the  following  thermo-elastic  formulation,  chemical  shrinkage  and  thermal  expan¬ 
sion  strains  were  taken  into  account  with  the  assumption  that  the  coeflUcients  of 
thermal  expansion  and  chemical  shrinkage  were  constant.  The  thermal  expansion 
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and  chemical  shrinkage  strain  were  dependent  on  the  change  in  temperature  and 
degree  of  cure,  respectively. 

The  equation  to  describe  the  instantaneous  resin  modulus,  specifically  in 
the  second  region,  was  a  rule  of  mixtures  model  based  on  the  degree  of  cure  [41] 
as  seen  in  the  following: 

■E'm  =  (1  -  oc)E^  +  aE^  +  7(1  -  «)«(£:“  -  (3.I) 

where  a  is  the  degree  of  cure  and  and  E^  are  the  fully  uncured  and  fully  cured 
matrix  moduli,  respectively.  7  is  a  fitting  parameter  corresponding  to  the  speed  at 
which  the  resin  inodulus  approaches  the  fully  cured  modulus.  The  resin  modulus 
approves  the  fully  cured  modulus  rapidly  if  the  value  of  gamma  is  high.  For 
the  following  analysis,  7  was  assumed  to  be  zero  to  represent  the  gradual  increase 
of  the  resin  modulus.  FVom  [41],  a  study  of  Equation  (3.1)  showed  reasonable 
agreement  with  experimental  results  with  7  =  0.  Similarly,  the  rule  of  mixtures 
was  used  to  compute  the  longitudinal  composite  stiffness,  Ei.  With  the  matrix 
modulus,  Ejn,  the  mstantaneous  composite  properties  were  computed  using  a 
self-consistent  field  micromechanics  model  from  Whitney  and  McCullough  [42] 
and  presented  in  Appendix  A. 
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3.4.1  Chemical  Shrinkage 

In  early  studies,  the  influence  of  chemical  shrinkage  on  residual  stresses 
was  neglected  by  assuming  that  the  chemical  shrinkage  occurs  in  the  process 
where  the  mechanical  properties  are  low  [1],  However,  recent  work  has  found  that 
the  chemical  shrinkage  has  significant  effects  on  the  residual  stress  development 
and  contributing  to  the  spring-forward  effect  in  many  composite  parts  [8].  The 
behavior  of  the  chemical  shrinkage  strain  during  the  curing  process  of  several 
thermosetting  resins,  including  AS4/3501-6  prepreg,  was  studied  by  Russell  (43]. 

He  reported  that  the  volumetric  shrinkage  can  be  modeled  as  a  linear  function  of 
degree  of  cure  in’ which: 


VkCa)  =  Kla  (3.2) 

where  is  the  total  volumetric  shrinkage  of  the  resin.  For  the  3501-6  epoxy 
resin  studied,  the  typical  value  of  is  5%  is  used  in  the  present  study.  Also,  the 
model  assumes  a  constant  coefficient  throughout  the  cycle.  The  relation  between 
the  volumetric  shrinkage  and  the  strain  due  to  shrinkage  is  seen  in  the  following 
expression: 

^ch  =  «m  =  1^1  +  14/i  —  1  (3.3) 

From  Equation  (3.3),  the  coefficient  of  chemical  shrinkage  («,„)  for  the  resin  can 
be  obtained  and  used  in  the  Whitney  and  McCullough  micromechanics  model  to 
determine  the  effective  coefficients  in  the  principal  directions,  («,<).  The  following 
expression  for  the  chemical  shrinkage  strain  is  obtained  based  on  the  change  in 
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degree  of  cure,  Aa{t). 

ef{t)=KjkAak{t)  (3.4) 

3.4.2  Thermal  Expansion 

Thermal  expansion  is  a  major  mechanism  contributing  to  the  residual  stress 
development  specifically  in  thick  section  composites  where  large  thermal  gradients 
exist.  In  addition,  the  thermal  expansion  coefficient  mismatch  between  the  resin 
and  the  fiber  produce  significant  stress  and  deformation.  The  expression  for  ther¬ 
mal  expansion  strains  is  modeled  as  a  linear  function  based  on  the  temperature 
difference.  The  coefficients  of  thermal  expansion  in  the  principal  directions  were 
assumed  to  be  constant  throughout  the  cure  cycle  and  obtained  using  the  mi¬ 
cromechanics  model  knowing  the  coefficients  of  the  constituents. 

efit)  =  <f>jkATkit)  (3.5) 

3.4.3  ABAQUS  Implementation 

ABAQUS  finite  element  software  was  used  for  the  stress  and  deformation 
analysis.  The  methodology  involves  the  utilization  of  the  temperature  and  de¬ 
gree  of  cure  distributions  obtained  from  the  Galerkin  FEM  cure  simulation  in 
conjunction  with  a  thermo-mechanical  model  in  ABAQUS.  In  order  to  model  the 
thermal  expansion  and  the  chemical  shrinkage,  a  user  subroutine,  UEXPAN,  was 
developed.  In  addition.  The  temperature  and  degree  of  cure  fields  were  obtained 
via  the  subroutines,  UTEMP  and  UFIELD,  respectively,  which  interpolated  the 


Table  3.1:  3501-6  Resin  Properties 
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K  (MPa) 

3.447 

BS  (MPa) 

3.447  X  10^ 

Kfcrtnfc  (%) 

5.0 

values  by  calling  and  reading  the  data  from  files  generated  by  the  cure  simula¬ 
tion.  The  user  subroutines  developed  for  ABAQUS  are  located  in  Appendix  B. 
To  model  the  degree  of  cure  dependent  modulus,  a  table  of  effective  composite 
properties  based  on  the  micromechanics  model  for  different  degrees  of  cure  was 
generated  in  Excel  and  used  in  the  ABAQUS  input  file. 

3.4.4  Nurnerical  Example 

The  residual  stress  and  deformation  development  of  the  curve-shaped  part 
simulated  in  Section  2.5.2  was  studied.  The  part  is  a  graphite/epoxy  composite 
with  unidirectional  fibers.  Figure  3.2  shows  the  mesh  detailing  the  coordinate 
system  used.  In  this  case,  the  fibers  are  oriented  along  the  x-direction.  The  values 
of  the  initial  and  final  modulus  of  the  resin  needed  to  calculate  the  instantaneous 
modulus  is  in  Table  3.1.  The  effective  composite  properties  were  determined 
using  the  instantaneous  cure  dependent  resin  modulus  and  the  constant  fiber 
properties  in  Table  3.2  in  a  micromechanics  model.  The  top  face  of  the  curve¬ 
shaped  part  was  constrmned  to  simulate  the  influence  of  the  mold.  A  uniform 
pressure,  167  kPa,  was  applied  to  the  bottom  face  of  the  part.  The  stre^  and 
strain  behavior  of  an  element  located  at  the  centermost  region  of  the  composite 
was  studied.  At  this  location,  the  element’s  temperature  goes  above  the  mold 


Table  3.2:  Elastic  Material  Properties  of  AS4  Fiber  and  3501-6  Epoxy  Resin 


Property 

AS4  Graphite 

3501-6  Epoxy 

-0.9 

57.6 

.^2 

7.2 

57.6 

Ku  (tie) 

0.0 

-0.0164 

K22 

0.0 

-0.0164 

^12 

0.2 

0.35 

1^13 

0.2 

0.35 

^23 

0.3 

0.35 

El  (GPa) 

206.8 

3.2 

E2  (GPa) 

20.68 

3.2 

Gi2  (GPa) 

27.58 

1.185 

Gn  (GPa) 

27.58 

1.185 

G23  (GPa) 

6.894 

1.185 
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temperature  because  of  the  heat  generated  from  the  cross-linidng  polymerization. 
To  study  the  mechanisms  driving  the  residual  stress  development,  the  following 
three  simulations  were  performed: 

(1)  Case  1:  Only  thermal  expansion 

(2)  Case  2:  Only  chemical  shrinkage 

(3)  Case  3:  Both  thermal  expansion  and  chemical  shrinkage 

Prom  the  simulation,  the  most  significant  stress  development  occurred  along  the 
fiber  direction.  The  maximum  normal  stress  in  the  x-direction  (an)  was  approx¬ 
imately  2  MPa  for  Case  3.  Figures  3.3  and  3.4  show  the  resulting  plots  of  the 
normal  stress  and  normal  strain,  respectively,  for  all  three  cases.  When  the  ther¬ 
mal  expansion  effects  are  only  taken  into  account,  the  strain  development  closely 
follows  the  temperature  profile  since  thermal  strain  is  driven  by  the  temperature 
difference.  Likewise,  the  chemical  shrinkage  strain  behavior  is  similar  to  the  de¬ 
gree  of  cure  profile.  The  stress  obtained  from  the  thermal  expansion  case  shows 
that  the  part  is  in  a  state  of  compression  reaching  a  maximum  stress  of -2.5  MPa. 
In  the  chemical  shrinkage  case,  the  stress  is  doubled  compared  to  Case  3  where 
the  effects  of  both  strains  are  taken  into  account.  The  stress  and  strain  profile 
for  Case  3  always  falls  in  between  Case  1  and  2.  Prom  this,  one  can  see  how 
the  mechanisms  of  chemical  shrinkage  and  thermal  expansion  compete  witli  one 
another.  As  a  result,  one  can  design  the  process  temperature  cycle  so  that  the 
thermal  and  chemical  strains  cancel  each  other  out  and  reduce  the  residual  stress. 
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The  shear  stress  and  strain  for  Case  3,  <^23  and  ^23  are  plotted  in  Figure  3.5.  The 
part  experiences  a  sharp  dip  in  the  shear  strain  (about  — 7.5xl0~®)  early  in  the 
cycle.  This  can  be  attributed  to  the  low  viscosity  of  the  resin  at  the  beginning 
of  the  cure  cycle.  The  contour  plots  of  the  displacement  and  Von  Mises  stress  at 
the  end  of  the  cure  cycle  is  presented  in  Figures  3.6.  A  maximum  displacement  of 
6.8xi0“6  m  was  experienced  on  the  bottom  face,  particularly  along  the  edges  of 
the  curvature.  In  contrast,  the  maximum  stress  was  on  the  top  face  of  the  part, 
mostly  in  the  areas  of  the  curvature.  Because  studies  have  shown  that  an  elastic 
analysis  overestimates  the  residual  stress,  a  fuUy  anisotropic  viscoelastic  model  is 


desired. 


STRESS  (Pa) 


TEMPERATURE 


STRESS  (Pa) 
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Figure  3.5:  Elastic  Shear  Stress  and  Strain 
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Figure  3.6:  Elastic:  Displacement  (a)  and  Von  Mises  Stress  (b)  Contours  at 
t=9960  sec 
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3.5  Viscoelastic  Formulation 

Prom  various  studies,  it  was  reported  that  an  elastic  formulation  overesti¬ 
mates  the  residual  stress  and  deformation.  Therefore,  the  viscoelastic  behavior  of 
the  resin  should  not  be  neglected  in  an  accurate  process  model.  The  model  used 
to  depict  the  viscoelastic  behavior  of  the  composite  was  based  on  the  work  by  Kim 
and  White,  [2],  [27],  [13],  and  [26].  The  general  expression  for  the  constitutive 
relationship  for  an  anisotropic  linear  viscoelastic  material  is  [2]: 

T,  X,  t,  H)  -  e^i{X,  t')df  (3.6) 

where  is  the  stress  tensor.  Ski  and  sfi  are  the  mechanical  and  non-mechanicaJ 
strains,  respectively  while  X  is  the  location  vector,  a  and  T  represent  the  degree 
of  cure  and  temperature  while  t  and  If  are  the  present  and  past  times,  respectively. 
Equation  (3.6)  can  be  expressed  in  matrix  form  to  yield  the  following: 

=  /oo  ~  ;■  =  1  -  6)  (3.7) 

Using  the  chemical  shrinkage  and  thermal  expansion  strain  analysis  in  Section 
3.4.1  and  3.4.2,  the  expression  for  the  non-mechanical  strains  can  be  described 
in  Equation  (3.8).  Typically  moisture  concentration  contributes  to  the  total  non¬ 
mechanical  strain,  but  for  this  study,  the  effects  are  ignored. 

ef{X,t)  =  <f>jk^Tk{X ,t)  -f-  Kjk^Cik{X ,t)  (3.8) 

Equation  (3.7)  can  be  further  simplified  into  the  integral  form  in  the  following 
equation  assuming  a  zero  stress  history  before  i  =  0  and  thermorheologically 


simple  behavior  at  a  constant  degree  of  cure. 


(Ti{t)  -  ^  Qij{a,T,X,^-  _  ef{X,l!)de  {i,j  =  1-6)  (3.9) 

where  the  reduced  times,  ^  and  are  described  as  the  following  with  07*  as  a 
degree  of  cure  and  temperature  dependent  shift  factor. 


«=/ 

Jo 


df 


t - - - df 

^  Jo  ar(a,T(fO) 


(3.10) 


3.6.1  Orthotropic  Material 


Prior  to  developing  the  model  for  the  stress  relaxation  modulus,  a  brief  sum¬ 
mary  of  the  viscoelastic  stress-strain  relation  for  orthotropic  materials  is  presented 
in  matrix  form. 
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where  the  components  of  the  Q  matrix  are  functions  of  the  relaxation  and  shear 
moduli  and  Poisson’s  ratios  corresponding  to  the  principal  directions. 
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Qn  =  (3.12) 

n  -n  —*'31  +  i^2i^'23  p 

Vl2  —  V21  —  - - -C'l 

Q,3=Q3.  =  ^fJ‘^^E. 

1  -  t/i3l/3i 

■*  '  '  ^^22  -^2 

^  ^  _  J^32  +  I^12l'31  n, 

V23  --  ^32  =  - 

^  _  1  -  V12V21  ^ 

V33  ““  - ^ - ^3 

Q44  =  <^23 
Q55  =  <J31 
Qee  =  Gi2 

A  =  1  —  t'12^'21  ~  ^*23*^32  ~  21/21  ^'32t'l3 

For  a  transversely  isotropic  material  in  the  y-z  plane,  Equation  (3.11)  re¬ 
duces  to  a  function  of  five  independent  time-dependent  engineering  properties, 
El,  E2  =  Ez,  V12  =  Viz,  i'23)  and  Gn  =  Gis.  The  reduced  stress-strain  relation  is 
seen  in  the  following: 
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3.5.2  Epoxy  Resin  Material  Characterization 

The  challenge  in  modeling  the  viscoelastic  behavior  is  describing  the  devel¬ 
opment  of  the  stiffiiess  matrix  based  on  the  processing  parameters.  Since  all  resins 
have  different  stress  relaxation  behavior,  proper  stress  relaxation  tests  must  be 
performed  for  each  resin  system.  Since  the  implementation  of  a  viscoelastic  resid¬ 
ual  stress  and  deformation  analysis  is  the  focus  of  this  study,  the  model  described 
in  [2]  for  the  stress  relaxation  behavior  of  3501-6  epoxy  was  used.  Kim  modeled 
the  stress  relaxation  modulus  as  a  generalized  Maxwell  model  with  a  discrete  ex¬ 
ponential  Prony  series.  The  general  form  of  the  cure  dependent  expression  is  as 
follows: 

Emipi, 0  =  E^{a)  +  E^(a)  Wu{a)exp  (3.14) 

where  is  the  fuUy  relaxed  modulus  and  E;;,  =  -  E^,  in  which  is  the 

unrelaxed  modulus.  Equation  (3.14)  is  simplified  by  assuming  15“,  E^,  and 
are  independent  of  degree  of  cure  and  temperature.  are  weight  factors,  are 

discrete  relaxation  times,  and  ^  is  the  reduced  time. 

=  0.15) 

W=1  J 

The  reduced  time  of  the  resin,  is  depicted  in  the  following  expression: 

In  Kim  and  White’s  experimental  studies,  the  slope  of  the  shift  function,  or, 
increased  as  the  degree  of  cure  decreased.  As  a  result,  the  following  linear  function 
of  temperature  with  degree  of  cure  dependent  coefficients  was  developed.  Figure 
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3.7  shows  a  plot  of  the  shift  factor  for  different  temperatures.  This  plot  displays 
the  increasing  slope  of  the  shift  factor  as  degree  of  cure  decreases. 

log{aT)  =  Ci(q!)T  +  02(0:)  (3-17) 

ci(a)  = -oiercp 

....  ^2(0:)  =  -T°ci(a) 

where  ai  =  IA/°C,  02  =  0.0712/°C',  and  T°  =  30^(7. 

To  determine  the  stress  relaxation  times,  the  trend  of  peaJc  relaxation  times 
for  different  degree  of  cures  was  studied.  The  peak  relaxation  time  is  short  at 
lower  degree  of  cure  and  lengthens  as  the  cure  advances.  Kim  observed  that 
the  peak  relaxation  time  increased  approximately  6  orders  of  magnitude  from  a 
degree  of  cure  of  0.57  to  0.98.  This  trend  was  then  compared  to  the  change  in  glass 
transition  temperature  (Tg)  where  temperature  is  in  terms  of  degree  Celsius.  Prom 
the  testing  of  samples  of  varying  degree  of  cures,  a  curve  fit  of  the  glass  transition 
temperature  based  on  degree  of  cure  was  developed  as  seen  in  the  following: 

Tg{a)  =  10.344  +  11.859a  +  178.04a2  (3.18) 

This  expression  was  normalized  with  the  glass  transition  temperature  of  a  ref¬ 
erence  degree  of  cure,  Oref  =  0.98,  to  obtain  the  chemical  hardening  function, 
/(«)■ 

=  /(a)  =  0.0536  +  0.0615a  +  0.9227a2 

Tg{aref) 


(3.19) 
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Table  3.3:  Master  Curve  Parameters  for  a®  =  0.98  Degree  of  Cure 


u 

Tu  (min) 

1 

2.922137x10^ 

0.0591334 

2 

2.921437x103 

0.0661225 

3 

1.82448x10® 

0.0826896 

4 

1.1031059x10’' 

0.112314 

5 

2.8305395x108 

0.154121 

6 

7.9432822x10® 

0.2618288 

7 

1.953424x10’’ 

0.1835594 

8 

3.3150756x10’2 

0.0486939 

9 

4.9174856x10’^ 

0.0252258 

With  the  assumption  that  the  glass  transition  temperature  as  a  function  of 
cure  behaves  similarly  to  the  peak  stress  relaxation  times,  the  relaxation  time  at 
different  degree  of  cures  can  be  expressed  as: 


log{T^{a)) 
log{TP{a  =  0.98) 


=  /(«) 


(3.20) 


Measuring  that  the  peak  stress  relaxation  time  for  a  =  0.98  is  10®  ®  minutes,  a  plot 
of  the  peak  relaxation  time  versus  degree  of  cure  can  be  seen  in  Figure  3.8.  With 


the  assumption  that  the  other  stress  relaxation  times  followed  the  same  behavior 


as  the  peak  stress  time,  Equation  3.20  was  rewritten  as: 


logMa))  =  fia)log{T^{aref))  (3.21) 

Using  the  models  of  the  shift  function  and  stress  relaxation  times  and  the  pa¬ 
rameters  of  the  reference  degree  of  cure,  are/  =  0.98,  a  semi-log  plot  of  the  resin 
modulus  [Eq.  (3.15]  versus  reduced  time  was  generated  as  seen  in  Figure  3.9.  The 
parameter  values  for  aref  are  located  in  Table  3.3.  The  values  used  for  and 
were  0.031  GPa  and  3.169  GPa,  respectively. 


3.5.3  Composite  Effective  Properties 


To  define  the  composite  material  properties,  the  relaxation  of  the  composite 
was  assumed  to  follow  the  same  behavior  as  the  resin  modulus.  Because  of  this, 

Equation  (3.14)  was  rewritten  in  the  following  equation  to  describe  the  composite 
modulus. 


T)  -  Q* (a, T)  +  Qij{a, T)  ^  Wu{a,  T)exp  [—^^7^]  (3.22) 

L  ra,(o)  J  '  ’ 

Again,  the  assumptions  that  the  weight  factors  and  the  composite  unrelaxed  and 
relaxed  moduli  are  independent  of  degree  of  cure  and  temperature.  Equation  (3.22) 


reduces  to: 


Qijioi,  T)  -  [— (3.23) 

a;=l  L  J 

where  is  the  fully  relaxed  modulus  and  =  Q“  -  in  which  Q?-  is  the 
unrelaxed  modulus.  is  the  elastic  stiffiiess  of  the  lamina  and  obtained  from 
the  elastic  micromechanics  model  in  Appendix  A.  The  value  of  the  frilly  relaxed 
stiffiiess,  ,  is  based  on  the  unrelaxed  stiffness  (QV  )  and  a  partition  factor,  Rj 
as  seen  in  the  following  expression: 


Q^  =  RfQ^  iO<Rf<l)  (3.24) 

Rf  describes  the  degree  to  which  a  material  completely  relaxes.  If  Rf  =  1,  then 
the  material  doesn  t  relax.  If  Rf  =  0,  then  the  material  completely  relaxes.  For 
the  graphite/epoxy  composite  studied  in  this  case,  i?/  =  7 

The  elastic  material  properties  of  the  graphite  fiber  and  epoxy  resin  used  in 
this  study  are  in  Table  3.2.  The  fully  unrelaxed  stiffnesses  and  other  composite 
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Table  3.4:  Elastic  Composite  Properties  of  AS4/3501-6 


Qn  (GPa) 

127.4 

Qi2  (GPa) 

3.88 

Q22  (GPa) 

10.0 

Q44  (GPa) 

2.57 

^11 

0.5 

4>22  =  <f>33 

35.3 

Ku  {lie) 

-167 

«22  =  «33  {lJ£) 

-8810 

properties  were  calculated  the  micromechanics  model  found  in  Appendix  A.  The 
results  are  presented  in  Table  3.4. 

3.6  ABAQUS  Implementation 

3.6.1  Isotropic  Viscoelasticity 

For  the  viscoelastic  analysis,  two  methods  were  attempted  to  simulate  the 
residual  stress  and  deformation  development  using  the  ABAQUS  finite  element 
software.  For  the  first  methodology,  the  viscoelastic  function  in  ABAQUS  was 
adapted  to  match  the  model  found  in  Equation  (3.15).  The  viscoelastic  function 
in  ABAQUS  is  limited  to  linear  isotropic  materials.  Therefore,  the  following  pro¬ 
cedure  can  only  be  used  to  model  the  resin  behavior.  This  procedure  also  required 
the  use  of  the  subroutines,  UEXPAN,  UFIELD,  and  UTEMP.  Another  method 
to  describe  the  viscoelasticity  of  an  orthotropic  composite  system  is  described  in 
Section  3-6.3.  In  ABAQUS.  the  basic  hereditary  integral  formulation  for  linear 
isotropic  viscoelasticity  is  the  following: 

cr{t)  =  f  2G{t  —  T')edf  + 1  f  K{t  —  T')^dt' 

Jo  Jo 


(3.25) 
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where  c  and  0  are  the  mechanical  deviatoric  and  volumetric  strains,  K  is  the  bulk 
modulus,  G  is  the  shear  modulus,  and  r  is  the  reduced  time.  The  reduced  time 
in  ABAQUS  is  expressed  as: 

(3.26) 

where  Ab  is  the  shift  function  modeled  as  the  Williams-Landell-Ferry  (WLF) 
equation.  Rather  than  using  the  WLF  equation  to  describe  the  shift  function, 
the  subroutine,  UTRS,  which  is  found  in  Appendix  B,  was  generated  to  replace 
the  WLF  equation  in  ABAQUS  with  the  shift  function  model  presented  in  Section 
3.5.2.  The  relaxation  functions  K{t)  and  G(t)  can  be  defined  as  an  exponential 
Prony  series  similar  to  Kim’s  viscoelastic  model  in  Equation  (3.15). 

=  Koo  +  ^  Kicxp  (3.27) 

G{t)  =  Goo  +  £  GiCXp  — 1  (3.28) 

In  Kim’s  model,  the  relaxation  times  are  functions  of  degree  of  cure.  The 
relaxation  time  in  ABAQUS  does  not  permit  field  dependent  relaxation  times.  As 
a  result,  a  time-temperature-field  superposition  scheme  had  to  be  developed  to 
shift  onto  the  master  curve  at  the  reference  degree  of  cure,  oiref  =  0.98.  In  order 
to  shift  onto  the  master  curve,  the  modulus  at  the  current  degree  of  cure  must 
equal  the  modulus  at  the  reference  degree  of  cure.  For  this  to  occur,  the  following 
expression  must  be  true: 


■®in(Q()0  —  EmiOrefy^') 


(3.29) 
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where  the  superscript  (’)  denotes  the  values  at  the  reference  degree  of  cure.  In 
Equation  (3.15),  the  only  quantities  that  are  cure  dependent  are  the  reduced  (^) 
and  relaxation  (tu)  times.  Therefore,  Equation  (3.30)  must  be  true  in  order  for 


Equation  (3.29)  to  be  satisfied. 


'rU0‘ref) 


(3.30) 


Using  Equations  (3.17),  (3.19),  and  (3.21)  into  Equation  (3.30),  log{aT{aTef))  is  a 


function  of  the  current  degree  of  cure  (a)  and  temperature  (T)  and  the  reference 


degree  of  cure  (a^e/)  as  expressed  in  the  following: 


log{Qq'(o(rgf))  —  /f  (O!,  Q!re/,  T)  (3.31) 

N  N 

H{a,  are/,  T)  =  Y^  logr^ia)  +  logoria,  T)  -  ^  logru,(ao) 

aj=l  cj=l 

A  plot  of  If  (a,  Ore/,  T)  versus  temperature  for  different  degree  of  cure  cases  is 
presented  in  Figure  3.10  to  illustrate  the  temperature-degree  of  cure  superposi¬ 
tion.  Equation  (3.31)  can  be  expanded  to  yield  the  following  noting  that  the 
temperature  dimensions  are  in  absolute  temperature,  K: 

(cjearp  - -  -  C2)(T' -  303)  = /f(a,  are/,T)  (3.32) 

Ore/  -  IJ 

Solving  for  the  temperature  needed  for  the  master  curve  shift,  T',  Equation  (3.32) 
becomes: 

T'  =  +  303  (3.33) 

Ciexp  -  C2 

As  a  result,  the  shift  factor  computed  in  the  user  subroutine,  UTRS,  and  supplied 
to  ABAQUS  is  the  following: 

ar(are/,T')  =  Aff(0(t'))  =  log~^  Ciexp  ( - - — -'j  -  C2  {T  -  303)  (3.34) 

LL  \o^ef-lJ  J 
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A  schematic  of  the  flow  of  variables  in  the  user  subroutine  to  obtain  the  proper 
shift  to  the  master  curve  is  presented  in  Figure  3.11. 


Figure  3.11:  Flow  of  Variables  to  Shift  to  Master  Curve  for  ABAQUS  Viscoelas¬ 
ticity 
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3.6.2  Numerical  Example 

For  this  model,  the  same  curve-shaped  part  as  seen  in  the  previous  exam¬ 
ples  was  used.  Since  the  model  can  only  be  used  to  describe  the  matrix  of  the 
composite,  the  simulation  was  conducted  in  order  to  gain  an  understanding  of  the 
behavior  of  a  viscoelastic  material.  Since  the  material  is  only  the  resin  matrix, 
many  of  the  mechanisms  giving  rise  to  residual  stress,  such  as  the  mismatch  in 
thermal  coefficients  and  mechanical  properties,  are  not  present.  As  a  result,  the 
stresses  evolved  were  only  on  the  order  of  0-135  kPa.  Despite  most  negligible 
stresses,  the  transverse  normal  stress  in  the  z-direction,  reached  a  peak  value 
of  135  kPa  as  seen -10  Figure  3.12.  At  the  beginning  of  the  cycle,  the  stress  profile 
is  dominated  by  the  thermal  expansion  strains.  As  the  composite  hardens,  the 
chemical  shrinkage  strains  start  to  elfect  the  stress  development.  The  behavior  of 
the  stress  profile  toward  the  end  of  the  simulation  is  similar  to  the  advancement 
of  degree  of  cure.  From  this  observation,  one  can  see  the  importance  of  chemi¬ 
cal  shrinkage  contribution  to  residual  stress  and  deformation  development.  The 
normal  strain  in  the  z-direction  is  presented  in  Figure  3.13. 

Again,  the  three  cases  conducted  for  the  thermo-elastic  analysis  were  per¬ 
formed.  Rrom  the  isotropic  viscoelasticity  simulation,  the  most  significant  stress 
development  did  not  occur  along  the  fiber  direction  as  previously  seen  in  the 
thermo-elastic  case.  Figure  3.14  shows  the  normal  stress  (an)  for  all  three  cases. 
In  the  case  of  thermal  expansion  only,  the  part  reaches  a  maximum  stress  of  -2.5 
MPa.  For  the  case  where  only  there  are  only  chemical  shrinkage  effects,  the  nor- 
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mal  stress  reaches  a  maximum  of  4.5  MPa.  The  significant  stress  development  in 
Case  1  and  2  shows  the  tremendous  infiuence  of  thermal  expansion  and  chemical 
shrinkage  on  the  residual  stress  and  deformation.  If  either  mechanism  is  ignored, 
an  incorrect  stress  state  will  result.  FVom  Figure  3.15,  the  strain  for  Case  3  is  very 
small  on  the  order  of  4x10-*.  The  shear  stress  and  strain  for  Case  3,  023  and  £23 
are  plotted  in  Figure  3.16.  The  sharp  dip  in  the  shear  strain  at  the  beginning  of 
the  cycle  as  seen  in  the  elastic  case  does  not  occur  in  this  situation.  The  contour 
plots  of  the  displacement  and  Von  Mises  stress  at  the  end  of  the  cure  cycle  is  pre¬ 
sented  in  Figures  3.17.  A  maximum  displacement  of  1.27x10-4  m,  which  is  higher 
than  the  displacement  in  the  elastic  case,  was  experienced  on  the  bottom  face, 
particularly  along  outer  edges.  In  contrast,  the  maximum  stress  was  on  the  top 
face  because  of  the  constraint  from  the  mold.  Since  this  example  only  modeled 
the  properties  of  the  resin,  the  stress  only  reached  a  maximum  Von  Mises  stress 


of  198  kPa. 
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Figure  3.12;  Isotropic  Normal  Stresses  in  Z-Direction  (<733) 
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Figure  3.13:  Isotropic  Normal  Strains  in  Z-Direction  (£33) 
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Figure  3.14:  Isotropic  Normal  Stresses  in  X-Direction(<Tii) 
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Figure  3.16:  Isotropic  Shear  Stress  and  Strain 
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Figure  3.17:  Isotropic:  Displacement  (a)  and  Von  Mises  Stress  (b)  Contours  at 
t=9960  sec 
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3.6.3  Orthotropic  Viscoelasticity 

With  the  limitations  of  the  viscoelastic  function  of  ABAQUS,  another  ap¬ 
proach  was  taken  to  describe  the  viscoelasticity  of  an  orthotropic  composite  sys¬ 
tem  as  modeled  in  Section  3.5.3.  A  user  subroutine,  UMAX,  (Appendix  B)  was 
developed  to  calculate  the  changing  stress  relaxation  moduli  at  each  time  step. 
For  this  subroutine,  the  Jacobian  matrix  of  the  constitutive  equation.  Equation 
(3.23),  is  computed  at  every  time  step  and  the  resulting  stresses  and  deformations 
are  provided  by  ABAQUS.  This  scheme  also  utilized  the  previous  subroutines, 
UTEMP,  UFIELD,  and  UEXPAN,  to  obtain  the  temperature  and  degree  of  cure 
history  and  to  include  the  chemical  shrinkage  and  thermal  expansion  strains. 

3.6.4  Numerical  Examples 

For  a  completely  orthotropic  viscoelastic  simulation,  the  user  subroutine, 
UMAX,  was  developed.  In  this  case,  the  residual  stress  and  deformation  were 
modeled  in  the  following  three  steps: 

(1)  In-mold  cure  phase  modeling 

(2)  Part  Cooling 

(3)  Mold  Constraint  Removal 

Like  the  isotropic  case,  the  significant  stress  development  occurred  in  the  Z- 
direction.  The  plots  of  the  normal  stress  and  strain  are  found  in  Figures  3.18 
and  3.19,  respectively.  The  maximum  normal  stress  attained  was  2  MPa,  which 
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is  similar  to  the  elastic  case  but  along  a  different  direction.  The  profiles  of  the 
isotropic  and  orthotropic  case  were  similar  along  the  Z-direction.  Prom  Figure 
3.19,  the  strain  along  the  Z-direction  is  dominated  by  the  thermal  expansion 
mechanism.  For  comparison,  the  plots  of  the  stress  and  strain  in  the  X-direction 
are  also  presented  in  Figures  3.20  and  3.21.  In  contrast  to  the  isotropic  case, 
there  was  some  stress  development,  along  the  fiber  direction  of  about  400  kPa. 
The  shear  stress  and  strain.  Figure  3.22,  were  the  same  as  the  isotropic  case. 
From  the  plots  of  the  contours  of  stress  and  displacement.  Figure  3.23  for  this 
case,  the  high  stress  portions  on  the  part  are  in  the  regions  where  the  part  starts 
to  curve.  The  ma^mum  Von  Mises  stress  was  23  MPa.  The  deformation  results 
show  that  there  are  large  deformations  along  the  outer  edge  of  the  face  despite 
the  fact  that  a  uniform  pressure  load  was  applied  to  it.  These  deformations  may 
lead  to  cracks  and  delamination. 

A  simulation  of  the  cool-down  and  mold  removal  was  conducted  for  this 
case.  Figure  3.24  shows  the  normal  stresses  up  to  the  end  of  the  cool-down  phase. 
When  the  cool-down  phase  starts  (at  about  14,000  sec),  the  normal  stress  in  the 
^direction  ramps  up  to  6.5  MPa.  This  shows  that  there  is  a  significant  stress 
development  during  the  cool-down  phase,  which  is  the  area  where  many  studies 
were  conducted.  However,  there  is  still  a  significant  stress  development  prior 
to  cool-down  that  should  not  be  neglected  in  a  simulation  of  the  manufacturing 


process. 
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Lastly,  the  mold  removal  was  simulated  by  only  constraining  four  nodes  on 
the  top  surface  of  the  geometry.  Figure  3.25  shows  the  deformed  shape  superim¬ 
posed  on  the  undeformed  geometry.  The  two  ends  of  the  geometry  spring-forward 
and  the  overall  thickness  of  the  part  has  reduced.  These  effects  lead  to  part  dimen¬ 
sional  instability.  These  results  convey  the  importance  of  accurately  simulating 
the  mahufacturing  process  of  composites  to  minimize  these  effects. 


STRESS  (Pa) 


-both 

-shrinkage 

-expansion 

“TEMPERATURE 


0  1000  2000  3000  4000  5000  6000  7000  8000  90 

TIME(aac) 

Figure  3.18:  Orthotropic  Normal  Stresses  in  Z-Direction  (^33) 
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Figure  3.24:  Orthotropic  Normal  Stresses  up  to  Cool-Down 
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Chapter  4 


RESULTS  AND  DISCUSSION 


A  three-dimensional  thermo-chemical  cure  simulation  for  polymer  matrix 
composites  is  developed  using  the  Galerkin  finite  element  method.  Several  nu¬ 
merical  examples  are  presented  depicting  the  spatial  and  temporal  gradients  of 
temperature  and  degree  of  cure  within  a  part.  The  developed  methodology  is 
capable  of  effectively  solving  large  mesh  problems  with  the  resources  of  a  per¬ 
sonal  computer.  The  results  obtained  from  the  present  methodology  were  in  good 

agreement  with  the  experimental  results  available  in  the  literature.  Temperature 

% 

and  degree  of  cure  gradients  are  illustrated  for  a  typical  composite  part.  The  mesh 
sizes  used  in  this  effort  are  relatively  small,  and  finer  time-discretization,  typically 
one  second,  is  used.  This  is  shown  to  capture  the  temperature  variations  due 
to  the  exothermic  reaction  in  progress  during  the  cure.  While  the  temperature 
fields  are  predictable  without  the  presence  of  internal  heat  of  reaction,  the  heat 
generated  typically  presents  a  control  problem  during  the  manufacturing  of  these 
components.  Prior  studies  indicate  the  use  of  simpler  physical  models  or  neural- 
networks  based  approximate  methods  to  capture  the  physics  of  the  problem  in  a 
predictive  control  strategy.  The  per  time-increment  computational  burden  is  in 
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the  order  of  200  milliseconds  for  one  second  of  simulation  time.  Therefore,  the 
three-dimensional  temperature  and  degree  of  cure  field  prediction  with  reasonable 
mesh  sizes  can  be  performed  and  embedded  within  the  process  control  without, 
perhaps,  resorting  to  approximate  methods.  The  in-core  memory  requirement  is 
quite  affordable  as  the  core  size  required  for  a  10,000  DOF  problem  is  estimated  to 
be  roughly  32  MB.  Thus  the  resources  of  a  typical  personal  computer  is  adequate 
to  perform  the  analysis. 

The  residual  stress  and  deformation  fields  induced  during  manufacturing 
are  driven  by  the  the  temperature  and  degree  of  cure  distributions.  Together 
with  an  appropriate  constitutive  relationship  describing  the  mechanical  behavior 
of  the  material  during  the  cure  process,  the  computed  temperature  and  degree  of 
cure  distributions  obtained  firom  the  cure  simulation  enable  the  prediction  of  the 
process  induced  residual  stresses. 

In  simulating  the  process-induced  stresses,  several  approaches  were  imple¬ 
mented  including  an  elastic  and  viscoelastic  formulation.  The  plots  of  stress  and 
strain  show  a  strong  dependence  on  the  temperature  and  degree  of  cure  history. 
These  stresses  can  lead  to  part  warpage  and  subsequent  failure.  Therefore,  in  mod¬ 
eling  the  manufacturing  process  of  composites,  one  must  not  neglect  the  process- 
induced  stresses  and  deformations.  In  addition,  the  curing  phase  showed  signifi¬ 
cant  stress  development.  As  a  result,  the  stress  development  prior  to  cool-down 
should  always  be  considered. 

Prom  the  preceding  study,  a  three-dimensional  model  and  numerical  simu- 
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lation  of  the  process-induced  stress  was  developed  based  on  a  constitutive  model 
from  the  literature.  The  simulation  is  adaptable  to  include  the  use  of  other  con¬ 
stitutive  relations,  material  models,  and  fiber  orientations  and  structures.  This 
open  simulation  structure  offers  great  flexibility  and  can  be  extended  to  simulate 
different  types  of  composites  and  manufacturing  processes. 

.  TJie  objective  of  the  study  was  to  develop  a  general  platform  to  model  and 
simulate  the  residual  stress  and  deformation  during  cure.  However,  since  the 
constitutive  model  used  was  obtained  from  the  literature,  a  more  comprehensive 
study  into  the  viscoelastic  behavior  of  composite  materials  is  necessary  in  order 
to  investigate  the  validity  and  reliability  of  the  results.  The  results  obtained  and 
presented  above  were  used  only  as  a  guide  in  the  study  of  the  mechanisms  driving 
the  residual  stress  development.  This  numerical  analysis  only  serves  as  a  platform 
to  study  the  constitutive  relationships.  Moreover,  the  material  behavior  differs 
from  material  to  material,  which  further  increases  the  complexity  of  the  problem. 
Extensive  experimental  testing  of  each  material  system  at  various  temperatures 
and  degrees  of  cure  is  necessary  to  accurately  depict  the  viscoelastic  behavior.  In 
OTder  to  fully  model  the  viscoelastic  behavior  of  composites,  experimental  testing 
of  the  composite  system  should  be  conducted  to  obtain  proper  curves  to  model 
the  shift  function  and  reduced  and  relaxation  times. 

Using  ABAQUS  finite  element  software  as  the  solver  for  the  stress  simulation 
had  several  limitations.  As  explained  earlier,  the  ABAQUS  viscoelastic  function 
could  only  be  implemented  for  an  isotropic  material  with  field  independent  relax- 
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ation  times.  As  a  result,  the  user  subroutine,  UMAX,  needed  to  be  developed. 
However,  employing  the  UMAX  subroutine  requires  extensive  knowledge  in  the 
behavior  of  the  material  and  the  software.  Xhe  inability  to  control  the  solving 
scheme  of  the  software  limits  the  range  of  problems  that  can  be  accurately  mod¬ 
eled.  Rather,  the  viscoelastic  behavior  should  be  simulated  with  a  self-developed 
code  to’  have  a  better  Control  over  the  solving  scheme.  With  a  special  purpose 
finite  element  code,  the  behavior  can  be  modeled  more  accurately  and  there  is 
more  reliability  in  the  results. 

Xhe  issues  involved  in  this  current  study  is  the  subject  of  future  work.  An 
investigation  intq  the  validation  and  reliability  of  the  results  will  be  conducted 
through  experimental  material  testing.  Experimental  testing  is  necessary  espe¬ 
cially  for  low  cure  states  (a  <  0.5).  Xhe  viscoelastic  model  used  from  the  litera¬ 
ture  only  had  experimental  validation  for  cure  states  above  0.5.  Also,  the  devel¬ 
opment  of  an  FEM  code  to  overcome  the  limitations  presented  by  ABAQUS  will 
also  be  investigated.  Modeling  other  manufacturing  processes,  fiber  orientations 
and  structures  is  also  the  subject  of  future  research.  In  addition,  the  viscoelastic 
study  will  be  applied  to  investigating  the  characterization  of  the  physical  aging 
and  degradation  of  polymer  matrix  composites  since  they  are  similarly  influenced 


by  time  and  temperature. 
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Appendix  A 


Micromechanics  Field  Model 


The  following  model  from  Whitney  and  McCullough  [42]  was  used  to  com¬ 
pute  the  effective  composite  material  properties  based  on  the  constituent  proper¬ 
ties  and  fiber  volume  fraction.  A  principal  material  coordinate  system  for  unidi¬ 
rectional  composites  as  seen  in  Figure  A.l  was  adopted  where  ”1”  or  ”L”  corre¬ 
sponds  to  the  direction  along  the  fiber  (longitudinal  direction).  The  ”2”  or  ”T” 
direction  is  the  transverse  direction  (perpendicular  to  the  longitudinal  direction). 
The  model  is  formulated  for  an  orthotropic  material  that  is  transversely  isotropic; 
therefore,  the  properties  in  the  ”2”  direction  are  the  same  as  those  in  the  ”3” 
direction.  The  thermal  expansion  and  chemical  shrinkage  coefficients  were  also 
computed  using  the  micromechanics  model.  The  subscripts  /  and  m  represent  the 
fiber  and  matrix  constituent  properties,  respectively,  and  Uf  is  the  fiber  volume 


fraction. 


Figure  A.l:  Coordinate  System  of  Composite  System  (Principal  Directions) 
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A.l  Engineering  Constants 

Isotropic  Plane  Strain  Bulk  Modulus  (Equation  used  to  calculate  the  plane 
strain  bulk  modulus  for  the  constituents-  hr/  and  krm  ) 


2(1  —  4t'j2-®'2 


(A.I) 


Longitudinal  Young’s  Modulus: 


=  El  =  EtfU,  +  Ei^(l  -  Vf)  -H  -  vLTffkrfkTmGTrmi'^  -  vj)vf 

{hr/  +  GTTm)kTm  +  {h^f  —  kTri^GrTm^f 


(A.2) 


Major  Poisson’s  Ratio: 


vlt  =  (*^12  =  1^13)  =  -^LTf){kTm  -  A:r/)Grrm(l  -  Vf)vs 

[{Erf  +  GTTm)kTm  +  {krf  -  kTm)GTTmJ^f] 


In-Plane  Shear  Modulus: 


glt = (g.2 = G.,) = 

[[^LT/  +  CrxTm)  (GiTf  “ 


(A.4) 


TVansverse  Shear  Modulus: 


GtT  =  (G23)  =  ^rrm[ferm(<77Tm  +  Gtt/)  +  ‘^GTr/GrTm  -b  krmjGTTf  -  GTrm)l'f] 

[h'm{GTTm  +  GrTf)  +  ^Grr/GrTm  —  {krm  +  2Grrm)(Grr/  -  Grrm)*'/] 


Plane  Strain  Bulk  Modulus: 


_  (kr/  +  GTTm)kTm  +  {Erf  -  kTTn)GTTmt'f] 
{krf  +  G-TTm)  -  {krf  -  kTm)vf 


TVansverse  Young’s  Modulus: 


Et  =  {E2  =  E3)  = 


-i - 1 - 1_  _L  ^LT 

4A:r  AGtt  El 
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Transverse  Poisson’s  Ratio: 


I/jT  =  (1/23)  = 


lEiikr  —  EjjEJx  —  Av^rpkxEjT 
^Exjkx 


A. 2  Expansional  Strains 


(A.8) 


The  forms  of  the  following  expressions  were  used  to  calculate  the  coefBcients 
of  thermal  expansion  and  chemical  shrinkage  with  (1^1  =  ef ),  etc. 

Longitudinal: 


-^1  -  ExfVf  +  Eirnil  “  Vf) 


(A.9) 


Transverse:. . 


E  _  E  _  r  I  /1  w^Lf^Lf^f  d"  ^I,tn'^in»(l  ^ 

=  €2  =  Kt-/-'/  +  -  ‘'/M— 


(A.10) 


i"' 


1 
.  i 


Appendix  B 


ABAQUS  User  Subroutines 


i 


1 

V 


■  ■? 


ooo  ooooaon  on  on  anno none 
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B.l  UFIELD/UTEMP 


C 

C 

c  Thermo-kinetic  viscoelastic  subroutines 
c 

/”  /  /  !  r  7/  in  / - \ 

/  /  /  /  _ /  /  /  /  _ 11/  /  /  /  / 

/  /_/  /  /  _//-.// _  /  / _  /  / — /  / 

\ _ ij  / _ // _ /  / _ // - / 


SUBROOTIhE  UFIELD  { FIELD ,  KFI  ELD ,  MSECPT ,  K3TE? ,  KIMC ,  TIME ,  NODE , 
1  COORDS, TEMP, DTEMF,NFIELD? 

INCLUDE  ‘  ABA_PARAM .  INC  ‘ 

DIMENSION  FIELD  {MSECPT,  NFISLD)  ,  TI14S(2),  COORDS  (3). 

1  TEMPaiSECPT)  ,  DTSI'^P  (MSECPT) 


C0^2‘I0N  /UTEMP_DATA/DATAVAL  (2,2,5000),  C_TIME  ( 2 )  ,  ATIM,  BTIM, 

1  AINCRTIM, IFTIME , ICODE , MNODE 
CHARACTER  FILENAJyi*256 
DATA  IFTIME/ 0/ 

CHARACTER  CONFPATH*256 , FILNAME*265 

do  similar  code  to  field  for  temperature,  except  read  temperature  from  the  v 
:  files. 

IF (IFTIME . EQ . 0 ) THEN 


READ  IN  THE  CONFIG  FILE 

CALL  GETENWAR(  *this_dir*,  CONFPATH,  LENV^/AR  ) 
FILNAME  =  CONFPATH  / / ’ \user . cfg * 

WRITE (*,*)*  CONFPATH ‘ , CONFPATH , FILNAME 
OPEJT  (  52  ,  FILE=FILNA-ME ,  5TATUS=  *  OLD '  ) 

READ  { 52  ,  *  )  ATIM,  BTIM 
READ { 52 , * ) FSTTIME , SECTIMS 
READ { 52 , * ) KNODS 
CLOSE (52) 

READ  IN  2  TIME  STEFS 


AINCRTIM  =  SECTIMS-FSTTIHE 
icode  =  1 

c 

c 

C_TIMS( icode)  =  FSTTIME 
ifinc  =  ATIM*FSTTIHS-3TIM 
call  make_f ilename ( f ilenam, ifinc) 
call  read_vtk(f ilenam) 
c 
c 
c 

ISINC  =  ATIM*SSCTIMB+BTIM 
C 

call  make_f ilename ( f ilenam, i sine) 
icode  =  2 

c_time{ icode)  =  SECTXME 
call  read_vtk (f ilenam) 


no  a  o  ooof'j  ooooon 
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IFTIME  =  1 
END  IF 
c 

c  See  if  the  next  step  has  to  be  read  in. 
c 

IF(TIME(2) .GE.C_TIME(icode) )THEN 
prev^tiine  =  C_TIME(icode) 

C_TIME  ( 1 )  =prev_tiine 
ifinc=ATIM*C_TIME{l)  +BTIK 
icode=l 

call  inake_filename(filenain,  ifinc) 
call  read_vtk(filenain) 

-if*(icode  .eq.  Dthen 
icode  =  2 
else 

icode  =  2 
end  if 

C_TIME< icode)  =  Prev^time  +  AINCRTIM 
ISINC  =  ATIM*C_TIME (ICODE) +BTIM 
call  nvake_filenaine(filenam,isinc) 
call  read_vtk(filenair.) 

END  IP 


c 

c  interpolate 
c 


c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 


IVRITE  (*,♦)'  EI^ERING  UFIELD  ’ 

WRITE  (*,*)  COORDS(l) ,COORDS(2) ,CCOKDS(3) 

SPRITE  ( *  ,  * )  •  WRITING  CURE  FROM  UFIELD  * 

AZ  =  (TIME{2)  -  c_tiine(l))/(c_time(2)-c_tiine(l)) 

FIELD (1,1)  =  dataval  (1, 2,NODE) (1-AZ)  +  dataval  (2, 2,NODE)  *  (A2) 

WRITE ( ♦ , * )  • KSTEP=  » , KSTEP ,  *  KINC= ' , KINC 

WRITE(**,*)  •Tn^2=:*,TIME(2)  ,  •  CTIME  ( 1)  =  *  ,  C  TIME  (1 )  , 

1  *CTIME(2)=* ,C_TIME(2) 

WRITE (*, *)  *AZ=' ,AZ 

WRITE  ( *  ,  ♦ )  •  DATAVAL1=:  *  ,  dataval  (1,2, MODE)  , 

1  *  DATAVAL2=:  *  ,  DATAVAL  (2,2,  NODE) 

WRITE(*, ♦)  *UFIELD' ,FIELD(1, 1) 


RETURN 

END 


/////-.  _//  _ //  I  /  /  _  \ 

/  /  /  /  /  /  /  _/  /  /L/  /  /_/  / 

//-./////  / _  /  /  /  /  / 

\ - /  /_/  / _ /  /./  /./„/ 


Load  in  the  temperature  variations  using  UTEM?. 


SUBROUTINE  UTEM? (TEMP , M3ECPT , KSTEP , KINC , Tim , NODE , COORDS ) 

INCLUDE  *  ABA^PARAM . INC ’ 

DIMENSION  TEMP  (MSECPT)  ,  TIME  ( 2 )  ,  COOPJ)S  ( 3 ) 

COMMON  /UTEKP^DATA /DATAVAL (2,2,5000), C_TIME { 2 ) , 

1  ATIM,BTIM,  AINCRTIM,  IFTIME,  ICODE,  I-INODE 

CHARACTER  FILENAM*256 

do  similar  code  to  field  for  temperature,  except  read  temperature  from  the 
tk  files. 

C 


non  o  o  o  o 


IF ( IFTIME . EQ . 0 ) THEN 


READ  IN  THE  CONFIG  FILS 

OPEI>r  ( 52 ,  F1LE=  •  C :  \RTM\Geu_ABQ\user .  cf  g  ’  ,  STATUS=  '  OLD  ’ ) 

READ ( 52 , * ) ATIM , BTIM 

READ ( 52 , ♦ ) FSTTIME , SECTIIffi 

READ{52,*)MNODE 

CLOSE (52) 

P^J)  IN  2  TIME  STEPS 


AINCRTIM  =  SECTIMS- FSTTIME 
iccde  I 
c 
c 
c 

c_TIME{icode)  =  FSTTIME 

ifinc  =  ATIM*FSTTIME+BTIM 

call  make^f ilename  (filenam,  ifinc) 

call  read_vtk (f ilenam) 
c 
c 

C  '*  i; 

IS INC  =  ATIH*SECTXME+BTIM 
c 

call  inake_f  ilename  ( f  ilenam,  isinc) 
icode  =  2 

c_tiine  ( icode)  =  SSCTIME 
call  read_vtk (f ilenam) 

IFTIME  =  1 
END  IF 
c 

c  See  if  the  next  step  has  to  be  read  in. 
c 

I?  (TIME  (2)  ,GT.C_TIKE  (iccde)  )THEI^T 
prev^time  C_TIMS { icode) 
c_t  ime  ( 1 )  =prev_t  ime 
if  inc=ATIM*C_TIME  (1 )  +BTIM 
icode^l 

call  inake„filenaine(filenajn,  ifinc) 
call  read_vtk(filenam) 
if  (icode  .eq.  Dthen 
iccde  =  2 
else 

icode  =  2 
end  if 

C_TIME (icode)  =  ?rev_time  +  AINCRTIM 
ISINC  =  ATIM*C_TIME  ( ICODE) -kBTIM 
call  inake_filenanie{filenain,isinc) 
call  read_vtk {f ilenam) 

END  IF 


c 

c 

c 

c 


c 

c 

c 

c 


interpolate 

V^JRITE(*,*)  ’VJRITING  TEl^PERATURS  FROM  UTEMP’ 

AS  =  (TIME(2)  -  c_time{l) )  /  (c_tinie{2) -c_time(l) ) 

TEMP(l)  =  dataval (1, 1,N0DE} ^ (1-AE)  +  dataval (2 , 1 , NODE) * ( AZ) 
WRITE  (*,*)  ’TIME2=:’  ,  TIME  (2)  ,  '  CTIME(1)  =  ’  .C__TIME(i)  , 

1  ’CTTME(2)  =  ’  ,C_TIME{2) 

WRITE  (*,  *)  'AZ=‘  ,-AZ 

WRITE ( * , * )  ' DATAVALl = ' , da taval (1,1, NODE ) , 


C  1  ’ DATAVAL2  =  * , DATAVAL (2,1, NODE ) 

C  VmiTE(*.^)  *TEMP* ,TEMP(1) 

RETURN 

END 


SUBROUTINE  inake_f  i  1  enetme  ( f  i  1  enarr.e ,  i  n t  e ) 

Character  filenajne*256 

CHARACTER  IFBUF*30 

WRITE {IPBUF,fmt=* (130) * ) INTE 

ILSN  ?  1  ,  . 

DO  I  =  1,  30 

IF(IFBUF{I:I)  .NE.  *  ')THEJT 
ILE^I  -  I 
GOTO  30 
ENDIF 
END  DO 

30  CONTINUE 

filename  =  •C:\RTM\Geu_ABQ\vtk 
l\rtm_'//ifbuf (ilen:30) //* ,vtk* 

WRITE (♦,*)*  filename : * , filename 

RETURN 

BTJD 

C  ' 

C 

C 

SUBROUTINE  read_vtk ( f ilenam) 

INCLX7DE  *ABA_PARAN.INC' 

COMMON  /UTEI«>_DATA/DATAVAL  (2,2,5000),  C_TIME  ( 2 )  , 
1  ATIM,  BTIM,  AINCRTIM,  IFTIME ,  ICCDE ,  lE^IODE 

character  f  ilenaxn*256,  BUF*256 

open(53, f ile^filenam, status=  *  old’ ) 

100  CONTINUE 

READ(53, fmt=* (A256) ’ , END=200 ) BUF 
i f (BUF ( 1 : 14 ) . EQ . ' SCALARS  TEMPI  ‘ ) THEN 
RSAD(53,fmt=* (A256) *)BUF 
IVRITE  (*,♦)»  READING  TEMPERATURES  ’  ,  f  i lenam 
DO  I  =  1,  MNODE 

READ ( 53 , ♦ , END=20  0 ) DATAVAL  f icode ,1,1) 

END  1X> 

END  IF 

if  (BUF(1:14)  .EQ.  ‘SCALARS  CUREI_  ‘  )  THEI4 
READ(53,fmt='  (A256)  *)BUF 
WRITE (*,*)*  READING  CURES ‘ , f i lenam 
DO  I  =  1,  MNODE 

READ ( 53 , * , £ND=20  0 ) DATAVAL (icode ,2,1) 

END  DO 
GOTO  300 
END  IF 
GOTO  100 
200  CONTINUE 

IVRITE(  ^,  ♦)  ’FILE  ENDED  ABRUPTLY’ 

CONTIN*UE 
CLOSE (53) 

RETURN 
END 


300 


no  o  o  n  o  o  oooooooo 
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B.2  UEXPAN 


fill  _ / 

////_/  I 

/  /_/  /  / _  / 

\ _ / _ //_  / 


1/  /  /  _  \/ 
///_//  /| 
/  _ / 


/ 

/ 

/  / 


!_/_/  /_/  !_/_/ 


/  / 
/ 

/ 


SUBROUTINE  UEXPAN(EXP.AN,  DEXPANDT, TEMP, TIME, DTIME,  PREDSF, 

1  DPRED, STATBV, CMNAME , NSTATV) 

INCLUDE  .■ABA_PARAM.INC‘ 

CHARACTER* 80  CMNAME 

DIMENSION  EXPAN ( * ) , DEXPANDT ( * ) , TEMP ( 2 ) , TIME (2 ) , PREDE? ( * ) , 
1  DPRED ( * ) , STATEV (NSTATV) 
coinoosite 

DATA  ALPll,ALP22,ALP33/0.5E-6,35.3E-6,35.3E-6/ 

DATA  ALPll , ALP22 , ALP33 /O . 0, 0 . 0, 0 . 0/ 


steel  thermal  expansion  coefficients 
DATA  ALPll,ALP22,ALP33/13.9E-6,13.9e-6,13.9E-6/ 

DATA  VSHRI?IK1  ,.V3HRINK2 , VSHRINXS /  -167e-6,-8810e-6,-8810e-6/ 
DATA  VSHRINKl , V3HRINK2 , VSHRINK3 /O. 0,0. 0,0.0/ 

DEXPANDT (1)  =  ALPll 
DEXPANDT  (2)  =  ALP22 
DEXPANDT (3)  =  ALP33 

EXPAN  { 1 )  =  DEXPANDT ( 1 ) *  TEMP ( 2 ) +  VSHRILTCl *  DPRED ( 1 ) 

SXPAN  ( 2 )  =  DEXPA.rroT  ( 2 )  *TEMP  ( 2 )  +  VSHRINK2  *  DPRED  ( 1 ) 

EXP.AN  { 3  )  =  DE:<PAl‘roT  ( 3  )  *TEMP  ( 2 )  +  VSKRINK3  *  DPRED  ( 1 ) 
i«rite{*,*)  ‘TEilP^TEMP 
c  write ( * , * ) . ' cure ' , predef i 1 ) , dpred { 1 ) 
c  writei*,*)  ' expan' , expand) , expan(2) , expan (3) 

RETURN 

ajD 


o  o  o  o  o  o 
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B.3  UTRS 


/  /  /  //_  _//  _  \  /  _ / 

/  /  /  /  /  /  /  /_/  /  \ _  \ 

/  /_/  /  /  /  /  _,  _/ _ /  / 

\ - /  /_/  /_/  |_|  / _ / 


SUBROUTINE  UTRS (SHIFT, TEMP, DTEMP,  TIME, DTIME, PREDEF, DPRED, 

1  STATEV , CMNAME , COORDS ) 

C 

INCLUDE  'ABA_PARAM.TNC‘ 

C 

CHAPACTER*80  a<!NAME 
C 

DIMENSION  SHIFT { 2 ) , TIME (2 ) , PREDEF ( * ) , DPRED ( * ) , STATEV ( * ) 

1  COORDS(*) 

C  user  coding  to  define  SHIFT(l)  and  SHIFT(2) 

DATA  TW1,TW2, TV/3, TV;4,TW5,’n';6,TW7,TW8,TW9/1.75E+03,1.75E+05, 

1  1 . 09E+07 , 6 . 62E+08 , 1 . 70E+10 , 4 . 77E+11, 1 . 17E+13 , 1 . 99E+14 , 

2  2.95E+16/ 

DATA  ALPREF/O.^S/ 

C  V/RITE  ( * ,  * )  •  ENTERING  UTRS  ' 

c  WRITE(*,*)  COORDS (1 ), COORDS (2) , COORDS (3) 

C  'WRITE  (*,*)'  TIME=  • ,  TIME  ( 1 )  ,  •  DTIME=  '  ,  IxriME 

TIME  { 2 )  =TI!<!E  ( 1 )  +DTIME 

C  WRITEC,')  •TIME=’.TIME(2),  •  DTIME= ', DTIME 
C  VmiTE{*,  *)  ’TEl-IPr:' ,TEMP,  •  DTEMP= ’  , DTEMP 
C  VmiTE{*,*)’Field=*,PREDEF(l) , •  DFIELD= •, DPRED (1 ) 

C  COMPUTE  CONSTANTS 

FALPP^F=0 .0535-^(0. 0615*ALPREF)  +  (C  .  9227*ALPREF*ALPREF) 

R1  =DLOG10  ( TV/I  *TW2  *TW3  •  TV/4  *TV;5  *TW6  *  TV/7  'TW8  ♦TV/9 ) 
R2=-1.4*DEX?(1.0d0/ (ALPRSF-l.OdO) ) -0.0712 
C  WRITE(*,*)R1 

C  WRITE(*,*)R2 

C  CALCULATE  T'  and  aT* 

FALPHA1=0.0536+ (0.0615*PREDEF(1) ) + (0.9227*PREDEF(1) ♦PREDEF(l) ) 
C  WRITE  (*,*)  FALPH.A1 

~  ALOGATl= { -1 . 4 ‘DEXP { 1 . OdO/ ( PREDEF { 1 ) -1 . OdO ) ) -0 . 0712 ) * (TEMP-303 ) 
C  (/RITE  (*,*)'  ALOGATl  • ,  ALOGATl 

TPRIME1=  {  (  (  (FALPHAl/FALPREF)  ♦?.!)  +ALOGAT1-R1 )  /R2 )  +303 
C  VmiTE  ( * ,  * )  •T?RII4E1=  •  ,  TPRIMEl 

SHIFT(1)=10*»  (R2’  (TP.»IMEl-303> ) 

PREDEF2 = PREDEF ( 1 ) +DPRED ( 1 ) 

TEMP2  =TEMP+DTEMP 
C  WRITE ( * , ♦ ) DTEMP 

C  V/RITS  (  • ,  * )  DPRED  ( 1 ) 

C  WRITE  ( * ,  ♦ )  PREDEF2 

FALPHA2  =  0 .0536+ (0 . 0615*PREDEF2)  +  (0 . 9227*PREDEF2  *FREDEF2 ) 
ALOG.AT2= (-1 .4*DEXP(1 . OdO/ (PREDEF2-1 . OdO) ) -0 . 0712) * (TEMP2-3C3 ) 

C  V/RITE  ( * ,  ♦ )  '  1 ' 

T?RIME2=  (  (  ( (F.\LPH.^2/FALPREF)  *R1 )  +ALOGAT2-R1)  /R2)  +303 
C  V/RITE  (*,*)  ■TPRIHE2' ,TPRIHE2 

SHIFT(2) =10** (R2* (TPRIME2-303) ) 

C  V/RITE  (  *  ,  * )  •  TEMP1=  ’  ,  TEMP 

C  WRITE (*,*)  •TEMP2=‘ ,TEMP2 

C  V/RITE (*,*)  •PREDEF1=’ , PREDEF (1) 
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WRITK(*,*)  •SHIFT{1)=' ,SHIFT(1) 
WRITE{*,*)  •SHIFT(2)=’ ,SHIFT(2) 
WRITE (*,*)' EXITING  UTRS ' 

RETURN 

END 


OO  O  O  O  oo  ooooooooooo 
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B,4 


UMAX 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


/  /  /  /  1/  // 

/  /  /  /  /L/  /  /I 

/  /^/  /  /  /  /  _ 

\ _ iJJ  1^ 


/-  _/ 
/  / 

/  / 

/_/ 


STRESS{NTENS)  -  PASSED  IN  AS  THE  STRESSES  TENSOR  AT  THE  BEGINNING 
OF  THE  INCREMENT,  SIGMA_I_J.  (THIS  SHOULD  BE  UPDATED  IN 
THIS  ROUTINE.) 


STATEV(NSTATV)  -  SOLUTION  DEPENDANT  STATE  VARIABLES. 

DDSPDE  CNTENS,  NTENS)  -  JACOBIAN  MATRIX  Of  THE  CONSTITUTIVE 


SSE  -  SPECIFIC  ELASTIC  STRAIN  ENERGY. 

SPD  -  PLASTIC  DISSIPATION. 

SCD  -  CREEP  DISSIPATION. 

RPL  -  •  VOLUMETRIC  HEAT  GENERATION  PER  UNIT  TIME  AT  'tHE  END 
OF  THE  INC. 

DDSDDT  (NTENS)  ^  -  VARIATION  OF  STRESS  INCREMENTS  WITH 
RESPECT  TO  TEMPERATURE. 

DRPLDE  (NTENS)  -  VARIATION  OF  RPL  WRT  STRAIN  INCREMENTS. 


DRPLDT  -  VARIATION  OF  RPL  WRT  THE  TEMPERATURE  INCREMENT. 

STRAN (NTENS)  -  TOTAL  STRAINS  AT  BEGINNING  OF  THE  INCREMENT. 
DSTRAN(NT.ENS)  -  STRAIN  INCREMENTS 

TIME  ( 1 )  -  VALUE  OF  STEP  TIME  AT  BEGINNING  OF  CXJRRENT  INCREMENT 
TIME  (2)  -  TOTAL  TIME  AT  BEGINNING  OF  CURRENT  INCREMENT 
DTIME  -  TIME  INCREMENT 

TEMP  -  TEMPERATURE  AT  START  OF  INCREMENT 
DTEMP  -  INCREMENT  OF  TEMPERATURE 

PREDEF  -  ARRAY  OF  INTERPOLATED  VALUES  OF  PREDEFINED  FIELD 
VARIABLES . 

DPRED  *  ARRAY  OF  INCREMENTS  OF  PREDEFINED  FIELD  VARIABLES. 

CMNAME  -  NAME  GIVEN  ON  *MATERIAL  OPTION,  LEFT  JUSTIFIED. 

NDI  -  NUMBER  OF  DIRECT  STRESS  COMPONENTS  AT  THIS  POINT. 

NSHR  -  NUMBER  OF  ENGINEERING  SHEAR  STRESS  COMPONENTS  AT 
THIS  POINT. 

NTENS  -  SIZE  OF  THE  STRESS  OR  STRAIN  COMPONENT  ARRAY 
(NDI  +  NSHR) 

NSTATV  -  NUMBER  OF  SOLUTION  DEPENDANT  STATE  VARIABLES. 

PROPS  (NPROPS)  -  ARRAY  OF  MATERIAL  CONSTANTS. 
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C  OTROPS  -  NUMBER  OF  MATERIAL  CONSTANTS. 

C  COORDS  -  ARRAY  CONTAINING  THE  COORDINATES  OF  THIS  POINT. 

C  DROT(3,3)  -  ROTAION  INCREMENT  MATRIX. 

C  PNEWDT  -  RATIO  OF  SUGGESTED  NEW  TIME  INCREMENT  TO  THE 

C  TIME  INCP.EMENT  BEING  USED. 

C  CELENT  -  CHARACTERISTIC  ELEMENT  LENGTH. 

C  DFGRD0{3,3)  -  .ARRAY  CONTAINING  DEFOPJyiATION  GRADIENT 

C  AT  THE  BEGINNING  OF  THE  INCREMENT 

C  DFGRD1(3,3)  -  ARRAY  CONTAINING  DEFORMATION  GPJUjIENT 

C  AT  THE  END  OF  THE  INCREMENT 

C  NOEL  -  ELEMENT  I.TRIBER 

C  NPT  -  INTEGRATION  POINT  NUMBER 

C  LAYER  -  L.AYER  NUI-IBER  (FOR  COMPOSITE  SHELLS 

C  AND  LAYERED  SOLIDS) . 

C  KSPT  -  SECTIQK  POINT  NUMBER  WITHIN  THE  CURP.ENT  LAYER 

C  K3TEP  -  STEP  NUMBE.R 

C  KINC  -  INCREMENT  HUMBER 


SUBROUTINE  UM.AT  (STRESS,  STATEV,  DDSDDE,  SSE,  SPD,  SCD, 

RPL , DDSDDT , DRPLDE , DRPLDT , 

STPAN,DSTEAN,TIME,  DTIMB,TEMP,DTEMP,  PREDEF,  DPRED,  CMN.AME, 
NDI ,  NSHR ,  rjTENS ,  NSTAT\,' ,  PROPS ,  NPROPS ,  COORDS ,  DROT ,  PNEWDT , 
CELENT , DFGRDO , DFGRDi , NOEL , NPT , LAYER , KSPT , KSTE  P , KINC ) 
c 

INCLUDE  •  A3A_P.ARAM .  INC  ’ 

C 

CHARACTER*  80  CMTLAKE 

DIMENSION  STP.ESS  (NTENS)  ,  STATEV (NSTATV)  , 

1  DDSDDE (NTENS,NTENS)  , DDSDDT (NTSNS)  ,  DRPLDE  (NTEJIS)  , 

2  ST.RAN(NTENS)  ,DSTRAN(NTENS)  , TIME  (2 ),  PREDEF (1)  ,DPRED(1)  , 

3  PROPS (NPROPS) , COORDS (3) , DROT (3, 3) , DFGRDO (3,3) , DFGRDI (3, 3) 

"^COMI^ON/TW/T'WREF  ( 9 )  ,  TW  ( 9 ) 

COmON/ WEIGHT /WEIGHT  ( 9 ) 

COMMON/STIFF/QI  JO  (6,6),  QIJINF  (6,6),  ILO.AD 
D.ATA  ILOAD/b/ 

COMMCH/SHIFT/ZETAdOOOO,  9.2)  ,K3,KI 
D.AT.A  KS,KI/0,0/ 

c  Nvrober  of  Elements  in  Part 

DATA  NELEM/448/ 

c  Reference  Degree  of  Cure 

DATA  ALPREF/0.98/ 

c  WRITE(*,*)  'ENTERING  UMAT’ 

I F ( ILOAD . EQ . 0 )  THEN 
DO  1=1,6 

DO  J=l,6 


”'1 

QI JO ( I , J) =PROPS ( 6 ♦ ( I-l ) +J) 

QIJINF(I,J)=PR0PS{36+6* (I-l)+J) 

END  DO 

END  DO  '  1 

DO  1=1,9 

TVffiSF ( I ) = PROPS {72+1) 

WEIGHT(I)=PROPS(81+I) 

END  DO  'i 

C  ZETA  INITITIALIZATICN 
DO  1=1, NOEL 

DO  0=1, NPT 
.  ZETA{I, J,1)=0.0 

ZETAd,  J,2)=0.0  i 

END  DO  ' 

END  DO 


IL0AE=1  1 

END  IF  ,  I 

C  COMPUTE  CONSTANTS  FOR  SHIFT  FUNCTION  ^ 

ALPHA= PREDEF ( 1 ) +  0 . 5  *  DPRED ( 1 )  ) 

TEMP_N=TEMP+0.5*DTEMP  ’ 

TIME_N=TIM3  ( 2 )  +  0 . 5  *  DTIME 

FAI,PREF=0 . 0536V (0 . 0615 •ALPREF) + (0 . 9227*ALPREF*ADPREF)  1 

FALPKA=0.0536+(0.0615*ALPHA)+(0.9227*ALPHA«ALPHA)  I 

c  write(*.+)  'alphas’ , ALPHA, •  TEMP_N= ' , TEMP_N  5 

C  wr i te ( * . • )  ' FALPRSF= ' , FALPF^F , ’ FALPHA= ' , FALPHA 


^  TO^l'^1^9*'^~'*  laultiplied  by  60  to  make  units  equivalent  to  seconds 

TW(I) =60* (10’* ( (FALFHA’DLOGIO (TWRE?(I) ) ) /FALPREF) ) 

END  DO 

c  write{*,*)  .TW(1) 

SHIFT=10*  * ( (-1.4  *DEXF (1 . OdO/ (ALPREF-1 . OdO) ) -0 . 0712 ) * (TEMP  N-303) ) 
C  WRITE (*,*)  SHIFT 

c  write(*,*)  dtime 

TF( (K3TEP.ME.KS) .OR. (KING .NE . KI) )  THEN 
DO  I=1,NSLEM 
DO  J=l,9 

ZETAd,  J.  l)=2ETAd,J,  1) +2ETAd ,  J,  2 ) 

END  DO 
END  DO 
KSsKSTEP 
K1=KINC 
END  IF 


ZETA (NOEL , NPT ,2)  =  il/SKIFTI*C. 5  *DTIME 
C  2ETA_NCW= (TIME_N) + ( (l/SKIFT) *DT1ME) 

ZETA_NOW=ZETA{NO?.L,N?T,  1)  +ZETA(tIOEL,N?T  2) 


TE.RN1=0 

c  write(*,*)  'weights' ,v/eiaht{l) 

c  DO  1=1,9 

C  TERK=WEIGHT  ( I )  ’DEXP  ( -ZETA_NO'.';/TV';  d )  ) 

c  write(*,*)  ' terms ', term 

C  TERM1=TEPJ11+TERM 

C  END  DO 

DO  1=1,6 

DO  J=l,6 


DDSD0E  ( I ,  J)  =QIJIMF  ( I ,  J)  +  ( QUO  ( I ,  J )  *TERI'I1 ) 

END  DO 
END  DO 

RPL=0 

DRPLDT=0 

DO  1=1,6 

STRESS (I) =0.0 
DDSDDT{I)=0 
DRPLDE(I)=0 
DO  J=l,6 

■  STRESS  (i)  =STRE3S  (I)  +  (DDSDDE'd,  J)  {DSTRAN(J)  )  ) 
EH'JO  DO 
EMD  DO 


c  if  ( (NOEL. EQ. 440) .AMD. (NPT.EQ.l) )  then 

C  WRITE  ( *  ,  * )  ’  tertqp  *  ,  KING ,  TEMP__N ,  TEMP 

c  'cure’,  KING,  PREDEFd),  DPRED(l) 

c  end  i f 

IF  { (NOEL. EQ. 572) .AMD. (NPT.EQ.l) )  THEN 
Vmi?E{*,*)  'TIMS=’ ,TIME(2) ,  DTIME 
;-^RITE  cemp  *  ,  KING ,  TEM?_N,  TEMP 

WRITS(*,*}  '’cure',  KING,  PREDEFd),  DPREDd) 
WRITE(*,679)  ('nf(I)  ,1=1,5) 

write 680)  (TO(I)  ,1=6,9) 

W’RITE  (’",*)  •  NPT=  ’  ,  MPT 

WRITE ( * , * )  *  KINC= '  , KING ,  ‘  ‘ ’ KSTE?= ’ , KSTEP 

WRITE  ( ^  ,  * )  '  Z;ETA_N0W=  ‘  ,  LETA^NOW 
WHITE{*, *)  *aT* , SHIFT 
write ( * , * ) ' ddsdde ’ 

W71ITE(*,677)  (  (DDSDDE  (I,  J)  .J=l,6)  ,1=1,6) 
fonnat*(‘A’,6G12.5/,  ’B’,6G12.5/,  ’C‘,6G12.5/,  'E*,6Gi2.5/ 
1  , »F' ,6G12.5/, *G' ,6G12.5/) 


V/HITE(*,*)  ’dscrNd)  =  ‘ ,  DSTRANd) 
v'‘;RITS{*,673)  (uSTRAN(I)  ,1  =  1,6) 

VMTH(*,*)  *strNd)  =  ‘ ,  STRAMd) 
v’iRITE(^,673)  (STRAN(I)  ,1=1,6) 

FORMAT(6(gl2.5/)  ) 

_  format ( ’TW’ , {5gl2 .5) ) 
format  ( ’T2  * ,  (4gi2 . 5) ) 

formate 'T1‘ ,G12.5/, ’T2' ,G12.5/, ’ T3 ‘ , G12 . 5/ , ’ T4 ‘ ,G12.5/ 

1  , 'TS* ,G12.5/, ’T6’ ,G12.5/, ‘T7' ,G12.5/, ‘TS’ ,G12.5/, ’T9‘ ,G12.5/) 
VJRTTEi*,^)  'stress  (1)  =  ’,  STRESS  (1) 

V/RITEt*  ,  573)  (STRESSd)  ,i=l,  5) 

END  IF 


c678 

c679 

c680 

c679 

c 

c 

c 


c 

c 

c 

c 

c 


c 


c 

c 

c677 

c 


RETURN 

EMD 


